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ABSTRACT 

The  analysis  of  the  Lloyd's  mirror  interference  of  multiple  narrowband  lines  from  a 
sound  source  obtained  with  a  single  sensor  is  considered.  A  straightforward  global 
least  squares  optimisation  is  discussed.  The  analysis  is  based  on  a  straight  line 
trajectory  model  of  the  source  motion,  and  is  configured  for  any  range  independent 
propagating  medium.  Two  methods  are  proposed,  each  being  based  on  a  discrete 
optimisation  on  a  predefined  parameter  grid.  These  are  a  simulated  annealing 
algorithm  and  a  guided  search  method.  The  least  squares  cost  function  is  defined  in  a 
manner  independent  of  the  intrinsic  shape  of  each  narrowband  emission.  Simulation 
data  based  on  an  iso-speed  sound  profile  are  used  to  explore  the  proposed  analysis. 
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Source  Localisation  using  Lloyd's  Mirror 
Fringes  on  Narrowband  Lines 


Executive  Summary 

This  report  considers  an  approach  whereby  the  interference  patterns  observed  in 
lofargram  (narrowband  frequency-time  plot)  displays  can  be  used  to  localise  the  range 
and  depth  of  an  underwater  source.  This  is  of  significance  since  the  approach  is 
suitable  in  the  analysis  of  separate  independent  lofargrams  such  as  obtained  from 
single  sonobuoy  sensors  where  otherwise  no  depth  information  can  be  obtained,  and 
the  range  information  is  potentially  highly  uncertain.  Having  localised  the  source,  the 
slant  range  can  be  used  along  with  the  measured  time  dependencies  of  the 
narrowband  lines  to  estimate  the  absolute  levels  of  the  source  over  an  extended 
azimuthal  range.  These  in  turn  are  required  when  estimating  detection  ranges  for  the 

source. 

Lloyd's  mirror  interference  patterns  are  often  seen  on  lofargram  plots  and  have  been 
considered  a  nuisance  when  interpreting  the  data.  The  structure  of  these  patterns 
however  is  intrinsically  dependent  on  the  range  and  depth  of  the  originating  source, 
and  in  principle  can  be  used  to  extract  these  quantities.  Whilst  the  source  broadband 
levels  are  often  too  low  to  allow  a  fully  developed  pattern  to  be  seen  and  analysed,  the 
significantly  higher  level  narrowband  lines  will  also  carry  the  interference  modulation. 
We  wish  to  extract  and  use  the  modulations  from  these  narrowband  lines. 

The  methodology  advanced  in  the  report  is  based  on  an  iso-speed  sound  profile, 
although  the  extension  to  more  realistic  profiles  is  straightforward.  We  assume  that  the 
source  moves  with  constant  velocity  and  depth  and  that  an  approximate  solution  is 
known.  Such  a  solution  would  typically  include  the  CPA  range  and  speed  of  the 
source  obtained  from  a  standard  analysis  of  the  narrowband  lines,  an  estimate  of  the 
operating  depth  of  the  source  and  the  sonobuoy  depth  setting.  This  defines  a  region  of 
parameter  space  which  along  with  a  cost  function  is  used  to  perform  a  multi¬ 
dimensional  least  squares  optimisation.  Because  of  the  complex  structure  of  the  cost 
function  we  propose  a  suitable  global  optimisation  strategy,  and  show  that  this  is 
capable  of  extracting  the  required  data.  The  end  result  is  the  extraction  of  the 
azimuthally  dependent  narrowband  source  levels. 
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1.  Introduction 


In  an  earlier  report  [1]  we  considered  how  realistic,  non-constant,  sound  speed  profiles 
would  affect  the  Lloyd's  mirror  interference  pattern  that  can  arise  with  moving  sub¬ 
surface  noise  sources.  Specifically,  we  were  interested  in  utilising  the  information 
inherent  in  the  interference  pattern  for  source  localisation  using  a  single  omni¬ 
directional  sensor.  Under  the  assumptions  that  the  source  is  at  a  fixed  depth  and  has  a 
known  constant  speed,  straight  line  motion,  the  resulting  time  dependent  interference 
pattern  allows  an  estimation  of  the  source  depth  and  range  from  the  sensor.  In  that 
study  the  structure  of  the  broadband  intensity  interference  pattern  evident  in  the  time- 
frequency  spectral  plot  (lofargram)  was  considered.  The  benefit  of  using  the 
broadband  interference  pattern  is  the  ease  with  which  the  frequency  dependence  of 
the  interference  fringe  structure,  which  provides  the  clues  to  a  unique  estimation  of  the 
source-receiver  relative  position,  can  be  incorporated  into  any  position  estimation. 

The  problem  often  arises,  however,  that  no  broadband  interference  pattern  is  visible  in 
the  lofargram.  This  may  result  from  the  source  broadband  emission  being  too  weak,  or 
from  a  weak  reflected  sound  wave.  As  narrowband  spectral  lines  can  be  significantly 
stronger  than  any  broadband  emissions,  they  may  still  exhibit  interference  effects,  and 
thus  be  used  to  estimate  the  source  position.  In  this  case  we  have  only  a  snapshot  of 
the  frequency  dependence  of  the  interference  fringes  at  a  few  select  frequencies. 
Although,  ideally,  data  obtained  at  one  frequency  is  sufficient  to  obtain  an  estimate  of 
the  source  relative  position,  the  inherent  signal  noise  means  that  at  least  two  well 
spaced  frequencies  are  required  if  the  position  is  to  be  determined  with  any  accuracy. 
This  report  considers  how  these  data  might  be  used  to  provide  such  an  estimation,  and 
as  a  consequence,  result  in  an  absolute  estimation  of  the  signal  strengths  and  beam 
emission  patterns  for  any  observed  narrowband  lines  arising  from  the  same  source.  It 
is  these  latter  quantities  that  are  of  ultimate  interest,  as  a  knowledge  of  them  is 
essential  in  determining  detection  probabilities  for  the  source,  or  for  source 
classification. 

The  aim  of  this  work  was  to  develop  a  tool  for  tertiary  analysis  of  acoustic  signatures, 
in  which  a  simple  computational  approach  to  the  analysis  of  narrowband  line  data 
could  provide  an  estimation  of  the  desired  positional  information.  It  might  be  possible 
to  develop  the  tool  into  an  operational  aid  for  estimating  the  position  of  an  acoustic 
source,  in  particular  the  depth  of  a  submarine.  The  present  paper  examines  the 
mathematical  basis  for  such  a  tool,  but  does  not  consider  the  operational  aspects. 
Realistic  simulated  data  were  used  to  test  the  tool  as  this  enabled  the  true  source 
parameters  to  be  known  and  did  not  require  this  paper,  which  is  essentially 
mathematical,  to  be  classified.  A  limitation  of  this  report  is  that  the  method  has  not  yet 
been  tested  on  real  data.  The  generally  poorer  quality  of  and  more  complex  structure 
in  observed  narrowband  lines  may  restrict  the  suitability  of  approaches  such  as  that 
proposed  here  to  extracting  any  source  localisation  information  intrinsically  available. 
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2.  Overview 

We  consider  the  situation  in  which  a  lofargram  containing  a  number  of  narrowband 
lines  is  available.  Ideally,  this  spectrum  should  be  calibrated  in  magnitude,  although 
this  is  not  essential  in  the  discussion  to  follow.  It  is  assumed  that  any  background 
contribution  has  previously  been  identified  and  calibrated  in  magnitude,  and  a 
(calibrated)  frequency  spectrum  is  available. 

The  source  trajectory  defines  a  range  of  aspects  or  "views"  of  the  source  as  seen  from 
the  receiving  position,  where  this  view  is  die  (sound)  signal  emission  strength  at  a 
particular  azimuth  and  elevation.  As  the  source  position  changes,  the  strength  of  the 
received  signal  will  vary  in  accordance  with  the  signal  strength  at  die  associated 
azimuth  and  elevation  angles,  and  with  any  spreading  or  propagation  losses  within 
the  medium. 

Over  some  time  period  about  the  time  of  the  closest  point  of  approach  (CPA),  each 
narrowband  line  will  have  some  temporal  line  profile  or  shape,  arising  from  the 
propagation  loss  and  the  source  signal  aspect  dependence.  Lloyd's  mirror  interference 
means  that  the  received  signal  strength  will  result  from  a  coherent  combination  of  the 
emission  strengths  at  two  separate  elevation  angle  values  (for  a  fixed  azimuthal  angle), 
along  with  any  reflecting  boundary  and  propagation  attenuation  effects.  Thus,  in 
general,  the  strength  of  the  received  signal  as  a  function  of  time  will  not  be  directiy 
related  to  the  source  signal  emission  at  any  particular  elevation  angle,  but  will  be  the 
result  of  signals  of  differing  intrinsic  strengths  coherentiy  interfering,  together  with 
any  attenuation  and  focusing  effects  of  the  propagation  medium  (as  controlled  by  the 
sound  speed  variation).  In  turn,  the  source  trajectory  determines  die  specific  mapping 
between  the  time  and  the  source  "view",  and  the  degree  to  which  the  propagation 
medium  will  contribute.  It  is  thus  likely  that  an  estimation  of  the  source  position  at 
CPA  (and  hence  for  straight  line  motion,  the  trajectory)  and  the  intrinsic  source  levels 
for  a  narrowband  signal  will  be  intractable  without  making  acceptable  simplifying 
assumptions,  and  being  able  to  identify  the  origin  of  at  least  some  of  the  observed 
signal  strength  variation.  It  is  on  this  latter  point  that  Lloyd's  mirror  interference  can 
be  of  benefit,  as  die  interference  modulation  of  the  narrowband  line  profile  is  often 
clearly  identifiable  by  the  approximately  regular  variation  in  the  signal.  The  temporal 
positions  of  the  minima  due  to  this  variation  is  strongly  dependent  on  the  geometry. 

In  the  analysis  to  follow,  a  number  of  simplifying  assumptions  are  made.  The  receiver 
depth  is  known.  The  source  is  taken  to  be  travelling  at  a  constant  velocity,  and  at  a 
constant  depth.  This  reduces  the  source  trajectory  description  to  three  parameters;  the 
horizontal  range  at  the  closest  point  of  approach  (CPA  range),  the  depth,  and  the 
speed.  The  medium  is  taken  to  be  homogeneous,  but  with  possibly  a  depth  dependent 
sound  speed.  The  reflection  as  described  by  a  single  reflection  coefficient  is  taken  to  be 
frequency  and  angle  independent.  This  assumption  is  reasonable  for  frequencies 
below  a  few  kilohertz,  at  least  for  surface  reflections  at  the  water-air  interface.  The 
source  sound  emission  (beam)  pattern  is  taken  to  have  only  an  azimuthal  dependence. 
In  the  method  to  be  discussed,  this  latter  point  primarily  is  only  relevant  for  the 
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interpretation  of  the  adjusted  received  signal  strength,  after  interference  and  range 
effects  have  been  removed,  although  a  strong  elevational  dependence  could  affect  the 
apparent  placement  of  the  interference  minima.  We  limit  ourselves  to  those 
narrowband  lines  with  a  sufficiently  high  signal  to  noise  ratio.  By  "sufficient"  we 
mean  that  the  interference  modulation  is  clearly  identifiable  above  die  random  signal 
level  fluctuations.  Generally,  this  will  be  accomplished  with  a  lower  signal  to  noise 
ratio  for  strongly  interfering  sound  waves  than  for  weakly  interfering  sound  waves 
since  the  interference  pattern  peak  to  trough  height  ratio  will  be  larger.  It  is  further 
assumed  that  we  have  prior  knowledge  of  the  source  speed,  the  time  of  closest 
approach  and  the  true  frequencies  of  the  narrowband  emissions.  Typically,  these 
would  be  obtained  from  a  Doppler  analysis  of  the  same  narrowband  line  data. 

The  work  discussed  below  is  presented  in  terms  of  surface  reflection  in  an  infinite 
depth  iso-speed  medium.  This  is  for  convenience  and  does  not  restrict  the  application 
to  either  bottom  bounce  and  non  spherical  spreading  applications.  In  principle  the 
modelled  loss  data  can  arise  from  an  analytical  expression  as  used  here,  or  from  direct 
numerical  evaluation.  No  attempt  has  been  made,  however,  to  assess  the  degree  to 
which  a  non-constant  sound  speed  profile  will  complicate  the  analysis.  This  is  an  area 
for  future  work. 

The  basic  approach  to  the  source  localisation  solution  is  to  recognise  that  the  Lloyd's 
mirror  interference  essentially  results  in  a  modulation  of  the  received  signal  line 
profile.  Considering  the  iso-speed  propagation  model,  the  received  intensity  for  two 
ray  interference  is  given  as 


V  *D  kr  '-'s  j 


hK 


R2r 


Rr 


,2nf . 


A 


l  +  Y2^  +  2Y^cos(^(RR-RD)+<t>0) 
Rr  Rr 


where 


a) 


where  Id  (Rd)  and  IR  (Rr)  are  the  direct  and  reflected  wave  intensities  (path  lengths),  f 
is  the  line  frequency,  Cs  is  the  (constant)  sound  speed,  r\  is  the  reflection  coefficient  and 
<J>o  is  the  phase  shift  arising  from  reflection  at  the  boundary.  Thus,  the  received  signal 
appears  as  approximately  the  expected  direct  wave  signal  intensity,  modified  by  a 
frequency  and  time  dependent  function  of  the  source  CPA  range,  depth  and  speed. 
Notice  that  the  azimuthal  and  elevational  dependence  are  subsumed  into  the  effective 
reflection  coefficient  y,  which  correctly  will  not  be  constant  in  time  unless  the  source 
signal  beam  pattern  is  isotropic  in  nature.  However,  because  we  have  explicitly 
removed  the  propagation  range  dependence,  it  will  remain  nearly  constant  over  some 
suitable  time  interval  about  the  CPA  time. 

Therefore,  for  any  number  of  possibly  harmonically  related  narrowband  lines 
originating  from  the  same  source,  the  received  signal  intensities  can  be  separated  into 
the  product  of  an  intrinsic  component  times  a  "common"  component.  The  latter 
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contains  the  information  concerning  the  source  position,  so  that  a  direct  estimation  of 
the  source  trajectory  based  on  any  number  of  narrowband  line  profiles  is  possible.  We 
indicate  the  "common"  component  as  such  because  it  is  frequency  dependent, 
however  this  dependence  is  calculable  once  a  model  of  propagation  is  assumed. 

The  problem  as  now  cast  is  to  identify  the  "common",  or  interference  modulation, 
pattern  in  the  narrowband  line  data  and  remove  the  intrinsic  component  of  each  line 
profile,  leaving  a  set  of  data  that  apart  from  the  line  frequency,  depends  on  the  same 
set  of  parameters.  A  global  search  parameter  estimation  can  then  be  performed  to 
provide  an  estimate  of  the  underlying  source  trajectory.  Without  any  further 
information  we  will  assume  that  the  data  are  normally  distributed  random  quantities, 
in  which  case  the  least  squares  estimator  provides  the  minimal  unbiased  estimate  of 
the  solution  point. 

As  indicated,  it  is  taken  that  the  source  speed,  time  of  CPA  and  line  frequencies  are 
known.  If  it  is  assumed  that  there  is  no  elevational  dependence  of  the  source  beam 
pattern,  the  reflection  coefficient  r\  will  be  constant  for  all  lines.  Thus,  a  global  search  in 
our  model  will  require  the  identification  of  three  parameters;  the  source  CPA  range 
and  depth,  and  the  reflection  coefficient.  The  small  number  of  parameters  makes 
feasible  a  grid  search  or  evaluation  over  a  defined  region  of  parameter  space.  In 
particular,  if  the  determination  of  the  reflection  coefficient  can  be  decoupled  from  the 
estimation  of  the  source  CPA  range  and  depth,  the  problem  is  well  suited  to  such  an 
approach.  Ideally,  this  is  the  case  as  the  source-receiver  geometry  defines  the  positions, 
and  the  reflection  coefficient  the  strength  of  the  interference  modulation  minima.  The 
introduction  of  noise  and  a  non-isotropic  beam  profile  will  alter  this.  If,  as  would  be 
most  likely,  the  source  speed,  CPA  time  and  line  frequencies  are  derived  from  a 
Doppler  frequency  shift  analysis,  an  estimate  of  the  CPA  range  will  also  be  obtained. 
The  accuracy  of  such  a  range  will  generally  not  be  high,  however  it  can  be  used  in 
determining  a  possible  search  grid  size  and  position.  The  cost  function  to  be  defined  in 
§2.3  effectively  decouples  the  reflection  coefficient  from  the  source  CPA  range  and 
depth  parameters,  so  we  consider  here  how  these  latter  two  quantities  may  be 
determined. 

As  the  aim  is  to  be  able  to  include  any  range  independent  propagation  loss  model,  we 
need  to  be  able  to  include  efficiently  general  propagation  loss  data  into  the 
calculations.  A  typical  propagation  loss  programme  would  provide  as  the  minimal 
output  tire  calculated  loss  as  a  function  of  horizontal  separation  distance,  for  an 
isotropically  radiating  source  of  unit  strength  (ie.  OdB  loss  at  1  m  separation)  at 
specified  source  and  receiver  depths.  Basing  any  parameter  search,  at  least  in  the  first 
instance,  on  a  grid  allows  such  data  to  be  directly  incorporated  in  a  reasonably 
efficient  manner.  For  a  given  receiver  depth,  at  each  source  depth  on  the  search  grid, 
the  loss  data  need  only  be  calculated  over  the  horizontal  ranges 

RMW=RCPA(mil1) 

Rmax  =  a/RcpaC™*)2  +(v-  At(max))2  ^ 
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where  V  is  the  source  speed  and  At  is  the  time  from  CPA,  to  encompass  all  the  grid 
CPA  range  values. 

This  suggests  that  the  calculations  proceed  in  three  steps:  a  region  of  valid  parameter 
space  is  defined  and  a  search  grid  is  constructed,  the  propagation  loss  data  is 
calculated  for  each  source  depth  point  on  the  grid  over  the  range  limits  set  in 
equation  2,  and  the  cost  function  is  evaluated  at  the  discrete  grid  points  such  that  the 
position  of  the  global  minimum  is  found. 

In  the  method  to  be  suggested,  this  approach  is  used  as  an  initial  estimation  of  the 
solution  point,  after  which  a  refined  value  is  obtained. 


2.1  Parameter  grid  determination 

A  typical  loss  curve  for  Lloyd's  mirror  interference  will  have  three  regions;  the  near 
field,  the  interference  field,  and  the  far  field.  In  the  interference  region,  the  loss  will 
oscillate  in  a  roughly  periodic  manner.  Both  the  source  depth  and  the  CPA  range  will 
effect  the  placement  of  the  interference  region.  As  either  of  these  parameters  is  varied, 
the  phase  difference  between  the  direct  and  reflected  waves  will  also  vary,  and  the  loss 
will  fluctuate  from  one  local  extremum  to  the  next  (eg.  from  a  local  maximum  to  a  local 
minimum)  The  step  sizes  then  must  be  selected  such  that  the  change  in  phase  remains 
less  than  n  radians.  Otherwise,  aliasing  of  the  sampled  loss  data  will  occur. 

The  variation  in  loss  derived  for  the  iso-speed  propagation  model  serves  as  the  basis 
for  determining  this  spacing.  From  the  equation  giving  the  received  signal  intensity  for 
two  ray  interference,  with  Id  =  Ir  ==  Io, 


I  =  If 


1  ,  T 
R2d  R 


2t] 


Rn  ’Ri 


-  cos(<J>  +  <J>0) 


27rf* 

with<j>  =  k(RR-RD),  k  =  — —  (3) 

L'S 


the  variation  in  loss  from  a  local  maximum  to  local  minimum  will  occur  over  a  range 
AR  such  that  Ac(>  =  n.  Expressing  the  ranges  in  terms  of  the  CPA  range  Ro  and  the  source 
and  receiver  depths  Ds  and  Dr  respectively,  we  have  for  surface  reflections 


RR=A/R2+(DS+DR)2+(V-At)2 

RD=VR2+(DS-DR)2+(V-At)2 


whence  we  obtain  the  partial  derivatives 
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1  _  5<t>  9Rr  5(j)  9RD  Ds+Dr  Ds-Dr 

k  9DS  9Rr  9Ds  9Rd  9Ds  Rr  Rd 
1  =  8Rr  t  9<|)  9RD  R0  R0  U 

k  9R0  9Rr  9R0  9Rd  9R0  Rr  Rd 

For  reflections  from  a  bottom  layer  at  a  depth  Zo,  (Dr+Ds)  is  replaced  with  2Zo-(Dr+Ds) 
in  equations  4  and  5.  Putting 


ADS  = 


AR0  = 


^9<j>v' 


\8DSJ 
'  9(j)  ^ 

<3R0> 


A<|> 

A(J) 


(6) 


wuth  A<|>  =  7T  gives  the  conditions  for  the  maximum  step  sizes  as 


ADS  < 


Cs  Ds  +  Dr 
2f^[  Rr 


Rd 


AR0<-^_ 

2fft 


(7) 


Recall  that  Rr  and  Rd  are  functions  of  the  two  parameters,  Ro  and  Ds.  The  variation  in 
the  maximum  depth  and  range  increments  is  shown  in  figure  1  for  Dr  =  100  m  and 
fo  =  150  Hz.  Equation  7  shows  that  the  required  grid  size  (number  of  elements) 
increases  as  the  square  of  frequency. 


In  producing  a  regular  parameter  grid  for  a  number  of  narrowband  lines,  we  need  to 
select,  for  the  maximum  line  frequency,  the  ADS  and  AR0  which  define  the  minimum 
along  the  surfaces  in  figure  1.  The  heavy  contour  lines  shown  therein  indicate  the 
appropriate  step  sizes.  Values  for  Ds  and  Ro  can  be  determined  iteratively  as  follows. 

Starting  at  the  smallest  source  depth,  a  sequence  of  source  depths  is  firstly  generated 
by  setting 


9k 

9R^ 


=  0 

=  Const 


9({> 

dD~s 


(8) 


and  solving  for  Ro.  This  wall  give  the  largest  k  (the  smallest  step  size)  for  the  specified 
fixed  Ds.  The  estimated  value  for  R0  is  then  used  to  produce  a  new  Ds  value  from 
equation  7a  and  the  procedure  is  repeated.  The  contour  of  the  minimum  ADS  and  a 
sequence  of  source  depths  is  thus  produced. 
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The  Ro  steps  are  estimated  at  the  calculated  sequence  of  source  depth  values  from 
equation  7b  by  calculating  the  increase  in  Ro  at  each  of  the  source  depths  and  selecting 
the  minimum  value.  This  leads  to  the  next  value  for  Ro  and  the  procedure  is  repeated. 
Thus,  a  regular  grid  built  from  the  minimum  step  sizes  is  generated. 

In  practice,  the  step  size  should  be  smaller  than  that  produced  from  equation  7.  The 
values  therein  give  the  maximum  sizes  to  prevent  aliasing  for  the  iso-speed 
propagation  model,  however  for  a  non-constant  sound  speed  profile,  the  loss  can  be 
expected  to  vary  over  a  smaller  range  [1].  Because  of  this,  the  only  sure  way  of 
ensuring  that  the  grid  spacing  is  sufficiently  small  is  to  actually  evaluate  the  loss  over 
the  grid  and  monitor  the  spacing  between  loss  maxima  and  minima.  The  grid  spacing 
must  be  small  enough  to  see  this  structure. 


2.2  Loss  data  determination 

According  to  the  above  scheme,  the  loss  data  are  evaluated  at  the  parameter  grid 
source  depth  values,  and  over  the  horizontal  ranges  identified  in  equation  2.  For  a 
given  CPA  range  (and  known  source  speed  and  CPA  time),  the  time  values  of  the 
experimental  narrowband  line  data  will  define  the  ranges  at  which  the  loss  data  must 
be  evaluated.  The  required  values  can  be  correctly  evaluated  from  tabulated  loss 
values  at  specific  ranges  as  long  as  the  tabulated  range  step  sizes  are  sufficiently  small 
to  allow  interpolation.  Thus,  we  choose  the  range  increments  for  each  given  Ds  value 
so  as  to  prevent  the  aliasing  of  the  loss  data.  We  follow  a  similar  method  to  that  for  the 
selection  of  the  Ro  values,  but  now  use 


and 


1  g(j)  _  R _ R_ 

k  SR  Rr  Rd 


where  R  =  ^R20  +(V- At)2 


-i 


to  evaluate  the  increment  in  the  horizontal  range,  R,  at  each  source  depth  value. 


(9) 


Again,  the  step  size  should  be  smaller  than  the  above  estimation.  When  sufficiently 
small  this  then  allows  linear  interpolation  of  the  calculated  loss  values  to  be  used  when 
estimating  the  loss  at  the  ranges  implied  by  the  experimental  times. 


2.3  Least  squares  cost  function 

The  critical  component  of  a  least  squares  optimisation  is  in  the  definition  of  the  cost 
function  to  be  minimised.  The  problem  in  the  current  situation  is  that  we  wish  to 
constrain  the  search  to  match  specific  features  of  the  experimental  line  profile  data,  viz. 
the  Lloyd's  mirror  component  of  the  signal  level.  In  general  it  can  be  expected  that  the 
received  signal  will  not  be  symmetric  about  the  CPA  time.  Further,  the  measured 
signal  will  always  have  an  attenuation  arising  from  the  spreading  loss.  It  must  be 
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assumed  that  there  is  no  a  priori  knowledge  of  the  shape  of  the  source  beam  pattern. 
The  feature  here  of  importance  is  the  positioning  of  the  Lloyd's  mirror  interference 
fringes.  Thus,  we  need  to  separate  the  signal  variation  into  that  arising  from  the  source 
beam  pattern,  ie.  the  attenuated  intrinsic  signal,  and  the  interference  modulating 
pattern.  Having  achieved  this,  the  parameter  estimation  can  be  performed  using  the 
modulating  function  alone. 

It  is  not  sufficient  to  directly  use  the  measured  line  profile  in  the  parameter  estimation. 
The  non-symmetry  and  spreading  attenuation  of  the  signal  will  introduce  a  bias  to  the 
parameter  estimation  as  the  sum  of  squares  will  be  sensitive  to  the  overall  shape  of  the 
line  profile. 

To  perform  a  separation  of  the  intrinsic  and  modulation  components,  it  is  observed 
that  typically  the  temporal  variation  of  the  interference  fringes  will  be  faster  than  the 
variation  of  the  intrinsic  signal  (excluding  the  random  noise  variation).  The 
interference  fringes  will  appear  as  a  pseudo  periodic  signal  in  many  cases.  A  simple 
rectangular  running  average  filter  of  an  appropriate  width  can  then  be  employed  to 
achieve  the  separation.  The  filter  window  size  must  be  chosen  according  to  the 
"periodicity"  of  the  interference  fringes,  and  will  differ  with  each  frequency.  To  handle 
the  edge  effects  of  such  a  filter,  the  filter  width  is  made  to  reduce  as  the  extremities  of 
the  data  sequence  are  approached.  In  these  regions  the  filtering  will  not  be  as  suitable, 
however,  owing  to  the  spreading  attenuation,  the  importance  of  these  regions  to  any 
sum  of  squares  function  can  be  reduced. 

The  most  straightforward  approach  to  determining  the  filter  width  at  each  frequency  is 
from  a  manual  evaluation  of  the  positions  of  the  interference  fringes.  In  the  process  to 
be  discussed  herein,  it  is  assumed  that  the  positions  of  the  interference  fringe  minima 
are  identified  prior  to  processing.  That  this  can  be  done  is  guaranteed  by  the  earlier 
assumption  that  the  interference  fringes  are  clearly  identifiable.  Knowing  the 
approximate  fringe  spacing,  a  filter  width  set  at  some  multiple  of  this  spacing  will  be 
adequate. 

Applying  such  a  filter  to  the  line  profile  data  will  generate  an  estimate  of  the  received 
intrinsic  line  profile.  This  estimate  will  improve  as  the  interference  modulation  better 
approximates  a  periodic  function.  Having  an  estimate  of  the  received  intrinsic  line 
profile,  the  interference  modulation  function  is  directly  extracted  from  the  data.  This 
modulation  function  will  be  independent  of  the  source  beam  pattern  and  any 
spreading  attenuation,  and  the  ratio  of  the  peak  to  trough  heights  will  be  directly 
dependent  on  the  degree  of  reflection. 

The  effect  of  applying  a  running  average  filter  to  various  interference  modulated  line 
profiles  is  shown  in  figure  2.  A  window  width  approximately  four  times  the 
interference  (temporal)  fringe  spacing  has  been  selected.  This  value  is  found  to 
adequately  reproduce  the  underlying  intrinsic  beam  profile.  Because  of  the  attenuation 
of  the  signal  with  the  time  from  the  CPA,  To,  smoothing  the  raw  signal  will  introduce  a 
bias  and  overestimate  the  signal  at  large  time  offsets.  This  will  become  more 
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pronounced  as  the  filter  width  is  increased.  Using  the  logarithm  of  the  raw  data  will 
reduce  the  variation  in  signal  strengths  over  the  filter  window,  and  thus  die  smoothing 
filter  will  lead  to  a  better  estimate  of  the  true  attenuated  signal  strength.  Figure  2a 
shows  the  effect  on  the  estimated  attenuated  intrinsic  line  profile  of  such  a  filter,  along 
with  the  true  line  profile  shape.  The  modulation  functions  are  shown  in  figure  2b.  This 
shows  that  where  the  interference  fringes  lack  a  strong  periodicity,  the  beam  profile 
estimation  is  inadequate.  The  representation  of  the  signal  at  times  close  to  To  depends 
on  the  phasing  between  the  direct  and  reflected  sound  waves.  For  destructive 
jj^^grfgj"gi\ce,  the  signal,  and  thus  the  modulation  function,  is  underestimated.  This 
should  have  little  effect  on  the  parameter  estimation  as  the  positions  of  the  fringes  in 
the  modulation  function  are  the  critical  factor  in  determining  the  source  position  and 
the  central  region  will  represent  only  a  small  part  of  the  entire  modulation  function. 
Likewise,  the  estimation  of  the  effective  reflection  coefficient  will  be  determined  from 
the  overall  variation  in  the  modulation  function,  and  thus  any  underestimation  should 
have  a  minor  effect  on  the  deduced  value. 

An  example  is  shown  in  figure  3  for  three  frequencies  for  which  the  phase  at  To  varies 
between  0  and  it.  An  isotropically  radiating  source  was  assumed,  with  a  CPA  range  of 
477  m,  source  and  receiver  depths  of  50  and  200  m  respectively,  and  a  source  speed  of 
5  m/s.  The  difference  between  the  actual  and  the  estimated  signal  profiles  depends 
somewhat  on  the  width  of  the  filtering  window  and  the  size  of  the  reflection 
coefficient,  t|,  but  for  destructive  interference,  varies  between  two  and  four  decibels  for 
the  above  geometry  where  0.1<n<1.0  and  the  filter  width  is  between  100  and  200 
seconds. 

Following  from  the  definition  of  equation  1,  we  define  the  interference  modulation 
function  for  a  frequency  U  as  the  time  dependent  loss  at  a  specific  line  frequency 
corrected  for  the  spreading  loss.  Thus,  for  spherical  spreading  we  have 

-n  9-Trf 

FM(i)  =  l  +  y2&  +  2y^cos(— ^(Rr-Rd)  +  <1»0)  (10) 

kr  kr  '-'s 

This  contains  the  (assumed  constant)  effective  reflection  coefficient  y.  The  least  squares 
cost  function  is  defined  in  terms  of  this  function,  and  the  received  signal  level  Ye  as 


with 


cf  =  Sw(i)(aE(i)-aT(i))2 

data 


W(i)  = ' 


YE(i) 


max 


K(i)} 


,  CJE(i)  = 


YE(i) 


YE(i) 


and  GT(i)  =  |FM(i)| 


(11) 


In  this  equation  the  double  bars  ( 1 G  | )  indicate  that  the  enclosed  quantity  is 
normalised  with 
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|G(i)|  = 


0(0 


for  function  G 


(12) 


This  is  done  as  only  the  shapes  of  the  experimental  and  modelled  line  profiles  are  of 

interest.  The  barred  ( G )  quantities  represent  the  filtered  line  profile  data.  Thus,  aE  is 
an  estimation  of  the  Lloyd's  mirror  modulation  function. 

To  include  numerically  evaluated  loss  data,  the  modulation  function  Fm  would  be 
replaced  by  the  calculated  loss  data,  but  this  still  requires  the  direct  ray  spreading  loss 
correction.  In  the  current  case  we  have  chosen  to  approximate  this  with  the  spherical 
spreading  correction  appropriate  for  an  iso-speed  medium.  A  model  dependent  form 
for  the  spreading  loss  experienced  by  a  single  direct  ray  could  be  evaluated 
numerically,  however,  correctly,  this  would  need  to  be  done  at  each  range  implied  by 
the  current  parameter  set. 

Each  datum  point  in  the  summation  of  equation  11  is  weighted  according  to  the 
relative  strength  of  the  estimated  received  intrinsic  signal  using  the  weighting  function 
W.  This  reduces  the  importance  of  the  lower  signal  to  noise  sections  of  the  data. 

When  multiple  lines  are  being  considered  concurrently,  each  narrowband  line  will 
generate  a  sum  of  squares  value  as  defined  above,  and  the  total  cost  function  is  defined 
as  the  geometrical  mean  of  the  individual  values,  giving 


C  = 


Vfreq 


for  N  frequencies 


(13) 


In  this  way  no  single  line  is  permitted  to  dominate  the  parameter  selection  unduly, 
and  a  fit  that  adequately  represents  all  the  data  is  obtained.  This  is  appropriate  as 
multiple  narrowband  lines  are  required  for  a  well  defined  solution. 


3.  Parameter  Estimation 

Without  further  information,  an  unconstrained  non-linear  least  squares  estimation  of 
the  parameters  is  necessary.  Figure  4  shows  the  cost  function  evaluated  for  a  single 
narrowband  line  over  a  region  of  parameter  space  stirrounding  the  known  solution 
point.  The  data  are  the  result  of  a  simulation  in  which  an  isotropically  radiating  source 
is  modelled.  The  simulation  geometry  has  a  source  depth  of  50  m,  a  CPA  slant  range  of 
500  m,  a  receiver  depth  of  200  m  and  a  reflection  coefficient  is  0.5.  The  source 
narrowband  line  frequency  is  150  Hz.  Evidently  the  cost  function  has  a  complex 
structure,  consisting  of  multiple  valleys  and  local  minima.  The  valleys  represent 
regions  of  approximately  constant  phase  difference  between  the  direct  and  reflected 
sound  waves.  The  correct  solution  point  is  identified  from  this  plot,  however  the 
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identification  of  this  point  through  a  function  minimisation  parameter  search  would  be 
critically  dependent  on  the  starting  point.  A  suitable  point  will  not  be  known  in 
advance.  With  a  frequency  independent  medium  the  length  of  features  on  the  cost 
surface  will  scale  inversely  with  frequency.  Thus,  the  higher  the  frequency,  the  closer 
the  initial  point  needs  to  be  to  the  solution  point  in  order  that  a  successful  parameter 
estimation  is  obtained. 

Two  approaches  may  be  taken  in  order  to  overcome  this.  A  statistical  search  of  the  cost 
function  minimum  can  be  performed  using  a  technique  such  as  a  simulated  annealing 
algorithm.  This  provides  a  constrained  random  walk  through  the  parameter  space, 
searching  for  the  global  minimum.  A  simulated  annealing  parameter  estimation 
requires  no  gradient  information,  however  the  estimation  can  become  expensive  in 
terms  of  the  number  function  evaluations  that  must  be  made.  Such  an  approach  can  be 
readily  applied  to  a  discrete  parameter  space  such  as  is  constructed  herein  and  an 
algorithm  is  discussed  further  in  appendix  A.  The  suggested  algorithm  is  found  to 
successfully  identify  the  global  minimum  for  the  simulation  cases  considered  in  this 
section.  Other  statistical  searching  methods  such  as  genetic  algorithms  also  may  be 
suitable,  but  have  not  been  considered. 

A  second  approach  is  a  guided  search  in  which  the  solution  region  of  the  parameter 
space  is  localised  by  the  direct  evaluation  of  the  cost  function  on  a  grid,  leading  to  die 
selection  of  an  appropriately  close  initial  point.  An  unconstrained  direct  non-linear 
parameter  search  can  then  be  performed  in  the  confidence  that  the  estimation  will 
converge  to  the  correct  solution  point. 

In  both  cases,  an  initial  region  of  parameter  space  must  be  defined.  A  Doppler  solution 
to  any  of  the  narrowband  line  frequency  shifts  will  produce  source  speed  and  CPA 
range  estimates.  Of  these,  the  range  value  will  usually  be  the  least  well  defined  as  it  is 
governed  by  the  rate  of  change  in  the  Doppler  frequency  with  time,  but  may  prove 
useful  in  providing  loose  bounds  on  the  CPA  range  solution  point.  No  information 
concerning  the  source  depth  can  be  gleaned  from  a  Doppler  analysis,  and  either  an 
estimate  based  on  other  information  or  a  scan  of  the  parameter  space  will  need  to  be 
made.  The  automatic  search  capacity  of  the  simulated  annealing  algorithm  allows  a 
generously  sized  parameter  space  to  be  assumed,  albeit  potentially  at  the  expense  of  a 
large  number  of  cost  function  calculations.  However,  with  the  present  methodology, 
these  calculations  still  require  the  pre-calculation  of  the  propagation  loss  data,  and  this 
may  be  the  limiting  factor  in  deciding  the  size  of  the  parameter  space  and  the 
resolution  of  the  search  grid.  Certainly,  with  a  large  parameter  grid  the  simulated 
annealing  optimisation  would  require  far  fewer  cost  function  evaluations  than  a  full 
evaluation  of  the  grid. 

For  the  guided  search,  and  without  further  information,  the  cost  function  will  need  to 
be  evaluated  over  an  arbitrarily  positioned  low  resolution  parameter  grid,  and  die 
regions  of  local  minimum  identified.  In  doing  this,  it  is  better,  following  the  reasons 
outlined  below,  to  consider  the  total  cost  function.  The  higher  frequency  data,  in 
particular,  are  important.  Because  the  number  of  grid  points  increases  as  the  square  of 
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the  highest  frequency,  even  with  a  coarse  grid  step  size  (eg.  near  to  the  limits  defined 
in  equation  7)  the  number  of  cost  function  evaluations  may  become  large  and  a 
compromise  between  the  size  of  the  parameter  space,  resolution  and  the  maximum 
frequency  will  have  to  be  made.  After  candidate  solution  points  have  been  located, 
these  may  require  a  closer  inspection  at  a  higher  resolution  so  as  to  ascertain  whether  a 
reasonable  solution  point  (the  CPA  range  and  source  depth)  is  exposed.  The  suitability 
of  a  candidate  solution  point  can  be  obtained  by  evaluating  the  modulation  functions 
for  each  of  the  narrowband  lines.  If  the  point  is  incorrect,  the  match  with  the  data  will 
become  increasingly  worse  as  the  line  frequency  differs  from  the  frequency  used  in 
deriving  the  estimation.  This  is  a  faster  method  than  evaluating  the  total  cost  function, 
which,  depending  on  the  range  of  frequencies,  may  require  a  large  number  of  grid 
points,  leading  to  a  heavy  computational  and  memory  demand. 

After  a  suitable  region  of  parameter  space  has  been  isolated  and  a  starting  point  is 
identified,  the  parameter  search  should  include  all  suitable  narrowband  line  data.  The 
search  does  not  need  to  be  constrained  to  discrete  parameter  values  at  this  point.  It  is 
important  for  two  reasons  to  use  the  total  cost  function  defined  in  equation  13.  The 
individual  modulation  functions,  and  hence  the  estimate  of  the  source  position,  will  be 
increasingly  more  sensitive  to  variations  in  the  source-receiver  geometry  as  the 
narrowband  line  frequency  increases.  Secondly,  the  total  cost  function  displays  a  more 
simplified  structure  than  do  the  individual  functions  defined  by  equation  11.  Using  the 
above  simulation  geometry,  figure  5  shows  the  total  cost  function  evaluated  for  six 
narrowband  lines  with  frequencies  between  50  and  750  Hz.  The  solution  point  is 
clearly  identified,  but  is  also  very  localised  in  parameter  space.  The  relative  flatness  of 
the  surface  in  comparison  to  figure  4  arises  because  apart  from  the  position  of  the 
minimum  at  the  true  solution  point,  the  positions  of  the  various  structures  on  the 
surface  are  frequency  dependent.  With  multiple  narrowband  lines,  the  variations  in 
the  individual  cost  functions  tend  to  cancel  except  at  the  true  position  where  all 
contribute  towards  a  minimum  value. 

It  should  be  emphasised  that  the  method  described  in  this  report  generally  would 
require  a  number  of  spaced  frequency  lines  to  generate  a  well  defined  solution.  As  the 
signal  to  noise  ratio  of  the  data  increases,  the  number  and  frequency  range  required  for 
such  a  solution  would  reduce  (eg.  with  ideal  data,  only  a  single  narrowband  line 
would  be  necessary  to  define  the  source  -  receiver  geometry).  No  assessment  of  the 
limitations  imposed  by  differing  signal  to  noise  ratios  has  been  made. 

In  the  above,  we  have  assumed  that  a  suitable  reflection  coefficient  is  being  used.  The 
next  section  discusses  how  the  dependence  of  the  cost  function  on  the  reflection 
coefficient  can  be  decoupled  from  that  for  the  source  depth  and  CPA  range.  Therefore, 
as  the  size  of  the  interference  fringes  is  directly  dependent  on  the  degree  of  reflection, 
we  can  estimate  a  value  directly  from  the  strengths  of  the  modulation  functions.  Using 
the  iso-speed  model,  the  bounds  on  the  strength  of  the  modulation  function  at  each 
frequency  can  be  calculated,  and  by  adjusting  the  reflection  coefficient,  these  can  be 
made  to  best  represent  the  strengths  of  the  modulation  functions  for  each  of  the 
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narrowband  lines.  An  example  is  shown  in  figure  6  for  some  of  the  same  simulation 
data  used  in  figure  5. 

3.1  Narrowband  analyses 

We  consider  two  Lloyd's  mirror  narrowband  line  analyses  of  simulated  data.  This 
allows  a  model  consistent  with  the  assumptions  outline  above  to  be  studied.  The 
source  travels  at  constant  speed  in  a  straight  line,  and  at  a  fixed  depth.  The  received 
signal  is  the  result  of  two  beam  interference  between  a  direct  and  a  surface  wave  with 
the  surface  boundary  having  a  constant  reflection  coefficient  of  0.5.  For  simplicity,  the 
propagating  medium  has  a  constant  sound  speed.  The  narrowband  line  signal  to  noise 
ratios  are  chosen  such  that  interference  pattern  is  clearly  identifiable. 

A  simplified  source  beam  pattern  with  only  an  azimuthal  dependence  is  modelled. 
Thus,  the  direct  and  reflected  beams  have  the  same  intrinsic  strength,  and  the 
estimated  source  beam  pattern  can  be  directly  compared  with  the  true  beam  pattern. 
The  estimated  beam  pattern  is  the  result  of  the  direct  division  of  the  measured  line 
profile  by  the  fitted  modulation  function  and,  correctly,  is  the  true  beam  pattern  in  this 
case.  We  take  die  azimuthal  dependence  to  be  a  sum  of  Legendre  polynomials, 

I0(M)  =  [ZaLPL(cose)j  (14) 


where  the  L=2  and  L=3  terms  represent  an  asymmetric  forward-aft  dipole  pattern.  The 
two  cases  to  be  considered  include  a  symmetric  and  asymmetric  forward-aft  pattern. 

The  source  and  receiver  depths  are  50  m  and  200  m  respectively.  The  slant  range  at 
CPA  is  500  m,  and  the  source  speed  is  4.9  m/s.  Six  narrowband  lines  are  considered;  at 
50,  150,  300,  400,  550  and  750  Hz.  Figure  7  shows  the  positions  of  the  broadband 
interference  minima  fringe  contours  and  these  narrowband  lines.  The  frequencies  are 
selected  so  as  to  include  a  differing  number  of  interference  minima.  All  source 
frequencies  have  the  same  emission  beam  pattern,  with  this  being  either  of  two 
choices.  The  symmetric  pattern  is  constructed  using  equation  14  and  the  coefficients 
(ai,  CC2,  (X3)  =  (1.0,  0.0,  0.75),  whilst  the  asymmetric  pattern  uses  the  coefficients  (ai,  a.2, 
a3)  =  (1.0,  0.5,  0.75).  The  emission  beam  patterns  for  these  coefficients  are  shown  in 
figure  8.  For  clarity  we  assume  that  the  source  does  not  emit  any  broadband  signals. 
Thus,  the  received  signal  comprises  the  tonals  and  a  broadband  background  noise. 

Using  the  above  parameters,  a  digitised  time  series  of  the  received  signal  is 
constructed  [2],  and  the  narrowband  data  are  extracted  using  the  programme 
AMADEUS  [3].  For  the  present  purposes  the  narrowband  data  in  the  form  of  received 
signal  strength  at  each  time  step  are  required. 
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The  received  signal  strengths  as  a  function  of  time,  and  with  the  background 
contribution  removed,  are  shown  in  figures  9a  and  9b  for  the  two  beam  emission 
patterns,  from  which  the  interference  fringes  are  clearly  identifiable.  The  purple  and 
red  curves  on  these  plots  respectively  indicate  the  estimated  and  true  received  intrinsic 
strengths  based  on  a  filter  width  of  four.  The  50  Hz  data  has  one  node  either  side  of  the 
time  of  CPA,  and  the  filter  produces  a  poor  estimate  of  the  true  shape.  The  150  Hz  data 
has  two  nodes,  and  the  filter  produces  a  result  which  is  marginal.  This  is  more  clearly 
shown  by  the  extracted  modulation  functions  in  figure  10.  The  lower  and  upper 
bounds,  evaluated  for  straight  line  propagation  and  a  reflection  coefficient  of  0.5,  are 
included  on  these  plots.  With  correct  filtering,  the  two  modulation  functions  should  be 
similar,  and  vary  between  these  limits. 

Equation  11  includes  an  optional  weighting  function  which  biases  the  cost  function 
towards  the  strongest  parts  of  the  received  narrowband  line  profiles.  The  effect  of  this 
factor  on  the  ability  of  the  cost  function  to  correctly  locate  the  solution  point  is  shown 
in  figure  11.  This  shows  the  individual  cost  functions,  evaluated  using  the  correct 
reflection  coefficient  of  0.5,  for  the  symmetric  and  asymmetric  beam  profile  150  Hz 
data  sets  respectively.  The  independence  of  the  cost  function  on  the  form  of  the  source 
beam  pattern  is  evident.  The  unweighted  cost  function  exhibits  a  significantly  more 
structure  than  the  weighted  function,  and  many  potential  solution  points.  Indeed,  in 
this  example  the  location  of  the  minimum  value  on  the  grid  is  considerably  removed 
from  the  correct  position.  When  the  weighting  is  applied,  the  structure  of  the  cost 
function  is  considerably  simplified,  and  the  correct  solution  point  is  identified.  The 
effect  on  the  total  cost  function  is  less  pronounced  because  of  the  averaging  done  in 
producing  this  value. 

Parameter  estimations  based  on  the  total  cost  function  evaluated  with  and  without 
data  weighting  indicate  that  the  weighted  estimations  produce  consistent  (and  correct) 
values  for  the  source  depth  and  CPA  range,  independent  of  the  value  of  the  reflection 
coefficient.  If  the  reflection  coefficient  is  also  varied,  the  weighted  estimation  also 
produces  the  correct  value  for  this  quantity.  Without  weighting  the  estimated  solution 
point  depends  on  the  chosen  value  for  the  reflection  coefficient.  Further,  when  this 
coefficient  is  allowed  to  vary,  values  dependent  on  the  initial  value  of  the  reflection 
coefficient  are  obtained. 

It  appears  that  the  weighting  of  the  data  is  useful  for  two  reasons.  The  structure  of  the 
individual  cost  functions  is  simplified,  which  improves  the  selectivity  of  the  starting 
point  for  any  search,  and  the  dependence  of  the  cost  function  on  the  source  depth  and 
range  at  CPA  is  better  decoupled  from  the  boundary  reflection  coefficient.  This 
separation  is  of  advantage  to  both  optimisation  methods  as  the  reflection  coefficient 
can  be  estimated  from  the  modulation  function  without  the  need  for  a  three 
dimensional  search.  In  particular,  it  is  useful  for  the  guided  search  where  a  good  initial 
point  must  be  defined.  Then,  a  nominal  value  for  the  reflection  coefficient  can  be 
assumed  (eg.  0.5  or  some  estimate  based  on  the  strength  of  the  modulation  function) 
whilst  a  suitable  initial  source  depth  and  range  are  determined.  A  concurrent 
optimisation  of  these  three  parameters  can  then  be  performed. 
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Parameter  correlations  are  shown  in  figures  16  and  17,  wherein  the  source  depth  and 
the  range  at  CPA  are  seen  to  be  strongly  correlated.  The  correlation  between  these 
parameters  and  the  reflection  coefficient  is  weaker  because  they  will  be  primarily 
sensitive  to  different  features  of  the  modulation  function  (viz.  the  fringe  positions  and 
the  fringe  peak  to  trough  variation  respectively).  The  weighting  of  the  data 
presumably  simplifies  the  cost  function  structure  and  better  separates  the  dependence 
on  the  parameters  because  it  biases  toward  the  data  about  the  time  of  CPA  where  the 
interference  modulation  will  be  clearest.  The  data  at  times  far  from  CPA  will  be  less 
affected  by  the  interference  modulation,  especially  after  the  last  interference  minimum, 
and,  with  the  signal  being  weaker,  will  be  more  susceptible  to  the  influence  of  noise. 

The  global  optimisation  fits  to  the  various  line  profiles  are  shown  in  figure  12.  Using 
these  model  curves,  estimates  of  the  source  emission  patterns  can  be  extracted,  and  are 
shown  in  figure  13.  The  emission  patterns  can  then,  under  the  assumption  of  a  purely 
azimuthal  dependence,  be  analysed  in  terms  of  the  three  term  Legendre  sum  of 
equation  14  to  produce  an  analytical  beam  pattern. 

The  simulated  annealing  algorithm  applied  to  these  same  simulations  provides  a  direct 
estimate  of  the  solution  point.  Considering  only  the  weighted  cost  function,  the 
algorithm  correctly  identifies  the  modelled  source  depth  and  CPA  range.  Three  runs 
were  made  in  which  surface  reflection  coefficients  of  0.25,  0.50  and  0.80  respectively 
were  assumed.  A  large  parameter  space  with  a  medium  resolution  grid  was  used.  The 
sampled  total  cost  function  surface  is  shown  in  figure  14  for  the  data  with  a  reflection 
coefficient  of  0.5.  The  sampled  cost  function  shows  how  the  annealing  algorithm 
concentrates  on  the  region  of  the  solution  point,  but  also  widely  samples  the  parameter 
space.  The  evolution  in  the  annealing  process  is  shown  in  figure  15.  The  number  of 
sampled  parameter  configurations  for  the  symmetrical  dipole  beam  pattern  was  1906 
with  the  annealing  requiring  13  cycles.  Following  the  discussion  in  appendix  A,  the 
hit  rate  was  2.1%  leading  to  1866  separate  function  evaluations,  or  2.2%  of  the  full 
parameter  grid.  Such  statistics  are  dependent  on  the  size  if  the  parameter  grid.  The 
present  size  was  315x274. 


4.  Uncertainty  analysis 

There  are  two  aspects  to  the  estimation  of  the  uncertainties  in  the  deduced  parameters. 
The  optimisation  is  performed  with  the  assumption  that  the  source  speed,  time  of  CPA 
and  receiver  depth  are  known.  The  stability  of  the  solution  point  to  these  needs  to  be 
established.  Of  these,  the  receiver  depth  and  time  of  CPA  should  be  well  known.  The 
speed  may  be  less  so,  with  this  depending  on  the  strength  of  the  narrowband  lines 
under  analysis.  Secondly,  an  estimate  of  the  uncertainty  in  the  least  squares 
parameters  arising  from  the  quality  of  the  fit  needs  to  be  determined. 

As  the  functional  dependence  of  the  location  of  the  minimum  in  the  cost  function  is 
unknown  in  general,  and  because  it  will  certainly  have  a  non-linear  dependence,  the 
easiest  way  of  assessing  the  dependence  of  the  solution  point  on  the  source  speed. 
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receiver  depth  and  the  time  of  CPA  (V,  Dr,  To)  is  with  a  numerical  variational  analysis. 
For  the  special  case  of  a  iso-speed  medium,  an  analytical  expression  relating  the 
dependence  can  be  derived.  The  iso-speed  case  does  not  produce  a  unique  solution 
point,  as  can  be  shown  from  a  consideration  of  the  expressions  for  the  phase, 
equation  3,  and  the  path  lengths,  equation  4.  Thus,  we  have  with  rearrangement. 


♦  =  k-(RR-RD) 


Rr  =  V(Ro+V2At2) 


1  + 


(Ds+Dr)2 
R2  +  V2  At2 
2 

7 


(15) 


which,  making  a  first  order  Taylor's  expansion  of  the  path  lengths,  gives 


<f>(t  =  T)=k(RR-RD) 

~  yk-^R2  +  V2AT2 


(Ds+Dr)2 

r2+v2at2 


(DS-Dr)21 

R2  +  V2AT2  J 


2kDsDR 

7Ro+V2AT2 


,AT sT-T0  and 


(Ds+Dr)2 

r2+v2at2 


«1 


Then,  when  the  phase  is  held  constant  at  <J>o  the  approximate  expression 


V  = 


R0  [~f2kDsDRV  1* 
AT  l  <|>0R0  J 


(16) 


(17) 


defines  the  relationship  between  the  source  speed,  depth  and  CPA  range,  and  die 
receiver  depth  for  a  constant  phase  at  each  time  step,  and  therefore  a  constant 
modulation  function.  A  non-unique  solution  point  is  also  expected  [1]  with  a  non¬ 
constant  sound  speed  profile. 

In  the  present  case  V,  Dr  and  To  can  be  assumed  to  be  independent,  and  the  scatter  in 
their  values  will  imply  a  scatter  in  the  solution  point.  The  independent  uncertainty 
bounds  on  the  source  depth,  range  at  CPA  and  reflection  coefficient  can  then  be 
determined.  Such  an  analysis  will  be  computationally  intensive  in  that  a  least  squares 
optimisation  starting  from  the  same  initial  point  needs  to  be  performed  for  each  triplet 
of  known  values,  and  the  set  of  triplets  needs  to  be  sufficiently  large  that  the  range  of 
each  quantity  is  well  sampled. 
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1  (385)  seconds  and  0.75  (4.9)  m/s  respectively.  The  solution  points  for  the  CPA  range, 
source  depth  and  reflection  coefficient  were  determined  for  a  sample  of  80 
configurations  drawn  from  these  distributions.  The  results  are  shown  in  figures  16 
and  17. 

The  expected  strong  positive  correlation  between  the  deduced  CPA  range  and  the 
source  depth  is  observed,  whilst  no  correlation  between  the  reflection  coefficient  and 
these  quantities  is  found.  The  magnitude  of  the  variations  in  the  source  depth  and 
CPA  range  is  consistent  with  equation  17.  The  independence  of  the  reflection 
coefficient  is  expected  as  the  phase  difference  between  the  direct  and  the  reflected 
sound  waves  is  independent  of  the  reflection  coefficient.  The  origin  of  the  strong 
correlation  is  the  source  speed.  The  scatter  in  the  deduced  source  depth  values  with 
the  variation  in  the  source  speed  is  presumed  to  arise  (equation  17)  because  of  the 
concurrent  independent  variation  of  the  receiver  depth.  A  particulary  definite 
relationship  between  the  source  speed  and  the  deduced  CPA  range  is  noted. 

The  resultant  CPA  ranges  and  source  depths  are  normally  distributed  about  the 
expected  mean  positions  (482  cf  477  m  and  50.7  cf  50  m),  with  2o  levels  of  76  and  6.9  m 
respectively.  Clearly,  the  least  squares  solution  point  is  sensitive  to  the  source  velocity. 
Any  accurate  estimation  of  the  source  trajectory  therefore  demands  that  this  quantity, 
in  particular,  is  well  defined,  as  equation  17  indicates  approximately  equal  fractional 
uncertainties  for  the  source  depth,  CPA  range  and  the  source  speed. 

Given  an  estimate  of  the  solution  point,  an  estimate  of  the  least  squares  uncertainty  in 
this  should  be  established.  Owing  to  the  non-linearity  of  the  problem,  only  an 
approximate  estimate  of  the  uncertainties  can  be  made.  The  cost  function  surface 
shows  the  correlation  between  the  source  depth  and  range  at  CPA,  for  all  other 
parameters  fixed.  A  contour  of  constant  value  (or  a  surface  when  three  parameters  are 
considered)  about  the  solution  point  indicates  an  uncertainty  bound  on  the  solution 
point.  The  confidence  level  to  be  assigned  to  any  contour  requires  a  knowledge  of  the 
probability  distribution  of  the  data.  Assuming  a  x2  distribution,  the  confidence  level 
can  be  calculated  from  the  F-distribution.  Thus,  the  contour  level  associated  with  a 
confidence  level  of  (1-a),  *£6,1  -  a) ,  is 


«*e,l-a)  =  C(0) 


F(p,n-p,l-a) 


(18) 


where  C(0)  is  the  magnitude  of  the  cost  function  at  the  solution  point  and  F(p,  n-p,  1- 
a)  is  the  F-distribution  for  p  variables,  (n-p)  degrees  of  freedom  and  a  confidence  level 
of  (1-a). 

An  alternative,  simpler,  approach  is  to  linearise  the  problem.  Then,  the  confidence 
levels  appropriate  for  a  linear  least  squares  optimisation  can  be  used.  In  summary,  the 
data  Y  are  assumed  to  be  represented  by  the  model  function  F(£,,0)  which,  using 
Taylor's  theorem,  is  expanded  about  the  solution  point  giving 
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with 


F(§u,e)*F;+£z?u-(e.-§.) 

i=l  ’ 

Fu°  =  F(^u,§)  and  Z°u=^|^ 


e  =  e. 


(19) 


In  equation  19  the  0  are  the  best  estimates  of  the  parameters  0,  F°  is  the  value  of  F 
evaluated  at  the  solution  point  and  the  Z°  are  the  linear  coefficients.  Each  of  these  is 
evaluated  at  the  data  points  £u.  Thus,  the  fitting  function  is 


(Y  -  Fu°)  =  £  Z°u  ■  (6j  -  6j )  +  eu 

i=l 


(20) 


The  linear  least  squares  uncertainty  estimate  is  obtained  from  (ZTWZ)_1  where  Z  is  the 
nxp  matrix  (n  data  pointsxp  parameters)  formed  from  the  Zj,u  and  W  is  a  diagonal 
weight  matrix.  Assuming  that  the  y2  distribution  is  appropriate,  the  independent 
uncertainties  80i  are  obtained  from  the  diagonal  elements  of  the  covariance  matrix 

Cv  =  (ZTWZ)  •  S2  •  t(n  -  p,  1  -  0.5  ct)2 

giving  80  =  Vdiag(Cv)  (21) 


where  S2  is  the  data  variance  and  the  t(n-p,l-0.5a)  is  the  Students  t-value  for  a 
confidence  level  of  (1-a). 

The  variance  is  estimated  from  a  sum  of  squares  (SS)  analysis,  using 

SSTOtal  =  SSF[X  +  SSPURE 

with  SST0TAL  =  XW(i)^|(i)  ,  SSm=Jw(i)  TE2(i)  and  SSPURE  =  S2-(n-p)  ^ 

i=l  i=l 

The  SStotal  is  the  total  variation  of  the  data  'F  about  zero,  the  SSfit  is  the  variation 
about  zero  due  to  the  fitted  function,  and  SSpure  is  the  variation  about  the  fitted 
function.  As  the  beam  profile  is  estimated  from  the  deviation  of  the  data  from  the 
modulation  function  (ie.,  apart  from  the  statistical  fluctuation,  the  fit  is  always  perfect), 
the  SSfit  is  estimated  from  the  smoothed  version  of  the  data,  'F . 

An  assessment  of  the  linearity  of  the  problem  about  the  solution  point  can  be  obtained 
from  a  comparison  of  the  uncertainty  estimates  derived  by  the  two  methods,  viz.  the 
contour  and  the  linearised  model.  Further,  for  a  linear  system,  the  cost  contour  given 
by  equation  18  should  be  symmetrically  placed  about  the  solution  point. 
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Uncertainty  estimates  using  these  two  methods  are  shown  in  figure  18  for  the  two 
parameter  (source  depth  and  CPA  range)  problem.  The  marked  box  indicates  the 
uncertainty  bounds  derived  from  e equation  21.  Both  tests  of  linearity  are  seen  to  fail. 
The  sizes  of  the  estimated  bounds  are  in  strong  disagreement,  indicating  a  breakdown 
of  at  least  one  of  the  assumptions  given  above.  The  appropriate  method  of  determining 
the  uncertainties  in  the  parameters  is  to  use  equation  18,  however  a  determination  of 
the  contour,  or  surface,  is  required. 

If  the  total  cost  function,  being  the  geometric  mean  of  the  several  separate  cost 
functions,  is  used,  the  above  analysis  is  incorrect  as  the  cost  function  is  no  longer  a 
pure  sum  of  squares.  As  this  will  be  the  function  to  be  minimised  in  the  algorithm 
described  in  this  report,  the  independent  linearised  bounds  can  be  used  as  quickly 
evaluated  guides  to  the  uncertainty  levels,  although  strictly  the  confidence  level  will 
not  be  known.  The  overall  uncertainty  in  the  variables  will  typically  be  dominated  by 
the  uncertainties  in  the  fixed  quantities,  in  particular  the  source  speed,  rather  than  the 
uncertainty  due  to  the  least  squares  optimisation,  and  some  sort  of  variational  analysis, 
such  as  described  above,  would  then  be  more  appropriate  in  defining  confidence 
levels. 


5.  Conclusions 

This  report  has  considered  how  to  utilise  the  Lloyd's  mirror  interference  pattern  often 
associated  with  measured  narrowband  signals;  in  this  case  those  signals  obtained  from 
a  single  omni-directional  sensor.  Rather  than  seeing  the  interference  as  a  nuisance,  the 
fringe  positions  can  be  used  to  estimate  the  source-receiver  relative  position; 
information  that  would  otherwise  be  unavailable  from  a  single  sensor.  The  present 
study  has  outlined  a  standard  least-squares  optimisation  approach  to  extracting  this 
information  from  a  set  of  narrowband  line  data.  Whilst  the  analysis  has  been  based  on 
an  iso-speed  medium,  the  structure  of  the  algorithm  is  equally  suitable  for  any  range 
independent  propagation  loss  model. 

We  propose  that  the  narrowband  received  signal  strengths  be  separated  into  two 
components  -  an  estimate  of  the  intrinsic  signal  strength  as  would  be  measured  in  the 
absence  of  interference,  and  the  interference  modulation  function.  The  interference 
modulation  is  directly  determined  by  the  source-receiver  geometry,  and  thus  this 
allows  multiple  line  data  sets  to  be  considered  simultaneously,  and  used  in  a  global 
least  squares  search  for  die  underlying  geometry.  Such  a  separation  is  properly  valid 
when  the  source  sound  beam  patterns  have  only  an  azimuthal  dependence,  and  the 
boundary  reflection  can  be  treated  as  glancing  angle  and  frequency  independent.  The 
least  squares  cost  function  is  derived  from  the  modelled  and  estimated  modulation 
functions.  This  modulation  function  is  estimated  by  passing  the  measured  signal 
strength  data  through  a  running  average  filter.  This  should  be  adequate  as  the 
interference  fringes  are  roughly  periodic,  at  least  for  those  data  near  to  the  time  of 
CPA.  For  those  narrowband  lines  where  the  interference  modulation  is  clearly  not 
periodic  (when  the  signal  frequency  allows  only  one  or  two  minima  either  side  of  the 


19 


DSTO-TR-0583 


time  of  CPA)  this  filtering  is  inappropriate,  and  the  method  cannot  be  used  on  these 
data. 

To  make  the  problem  tractable,  the  source  is  assumed  to  move  at  a  constant  velocity 
and  depth.  This  reduces  the  parameterised  trajectory  description  to  four  variables  -  the 
source  speed  and  depth,  the  range  from  the  receiver  at  closest  point  of  approach  and 
the  time  of  closest  approach.  We  assume  that  the  speed  and  time  of  CPA  are  known, 
along  with  the  fixed  receiver  depth.  Then  the  optimisation  is  required  to  determine  the 
source  depth  and  CPA  range  as  well  as  a  boundary  reflection  coefficient.  An  unique 
solution  cannot  be  determined  for  the  iso-speed  medium  in  that  the  source  speed,  the 
source  and  receiver  depths  and  the  CPA  range  are  coupled.  Thus,  the  accuracy  of  any 
estimated  solution  point  will  be  strongly  dependent  on  the  least  well  determined  of 
the  "known"  quantities.  This  will  typically  be  the  source  speed,  which  can  be  extracted 
from  a  Doppler  analysis  of  the  frequency  dependence  of  the  same  narrowband  line 
data. 

We  propose  that  the  optimisation  is  based  on  a  discretisation  of  the  source  depth  and 
CPA  range  parameter  space.  Ideally,  the  reflection  coefficient  is  uncoupled  from  these 
parameters,  although  signal  noise  and  elevation  dependencies  will  weaken  this.  Data 
weighting,  which  effectively  biases  towards  the  CPA  time  where  the  data  is  most 
clearly  affect  by  the  interference,  improves  this  separation.  The  grid  spacing  must  be 
based  on  the  anticipated  variation  of  the  propagation  loss  with  die  source  depth  and 
CPA  range,  and  the  iso-speed  model  allows  an  analytical  estimation  of  this.  With  a 
suitably  small  grid  cell  size,  no  further  useful  information  concerning  the  depth  and 
CPA  range  will  be  obtained.  The  structure  of  the  cost  function  on  such  a  grid  is  found 
to  be  complex  and  no  single  minimum  can  be  expected.  Rather,  a  series  of  valleys,  each 
representing  an  approximately  constant  phase  relationship  between  the  direct  and 
reflected  waves,  will  be  present.  Two  problems  are  posed  by  this.  Firstly,  the 
appropriate  bounds  on  the  parameter  space  must  be  ascertained  and,  secondly,  the 
manner  of  the  optimisation  needs  to  be  considered.  The  parameter  space  bounds  can 
realistically  only  be  decided  from  other  information,  such  as  the  estimated  CPA  range 
produce  in  a  Doppler  analysis,  or  from  practical  limits  on  the  source  depth. 

The  structure  of  the  cost  function  is  considerably  simplified  by  including  multiple 
narrowband  line  data  sets,  covering  a  wide  frequency  range.  The  total  cost  function  is 
taken  as  a  geometrical  mean  of  the  individual  line  data  cost  functions.  In  a  frequency 
independent  environment  such  as  we  have  here,  the  size  and  positions  of  the 
structures  in  the  cost  function  will  scale  with  frequency  (ie.  they  are  strictly  phase 
dependent),  and  therefore  any  combined  cost  function  will  tend  to  accentuate  the  true 
solution  point  minimum  whilst  averaging  out  other  variations  in  the  cost  function 
surface.  This  will  help  in  obtaining  the  proper  solution,  however  the  problem  of 
multiple  minima  is  not  removed.  Two  methods  are  suggested.  The  first  is  to  apply  a 
simulated  annealing  algorithm  to  the  optimisation.  Such  algorithms  are  particularly 
suited  to  detecting  global  minima  amongst  multiple  local  minima.  The  second 
approach  is  to  iteratively  search  the  parameter  space  looking  for  potential  solution 
points  and  testing  these  against  the  collection  of  the  line  data.  We  propose  that  this  is 
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done  using  the  total  cost  function,  and  with  a  coarse  grid  cell  size.  The  higher 
frequency  data  will  improve  the  sensitivity  to  the  parameter  values,  although  at  the 
expense  of  a  smaller  grid  cell  size,  and  potentially  a  longer  computation  time.  Indeed, 
if  the  bounds  on  the  parameters  is  not  large,  it  may  be  appropriate  to  evaluate  the  total 
cost  function  over  the  full  grid  and  directly  isolate  the  solution  point.  Once  a  suitable 
region  is  isolated,  a  direct  search  optimisation  algorithm  can  be  employed  to  provide 
the  appropriate  solution  point.  The  important  point  is  that  the  weighted  sum  of 
squares  cost  function  based  on  the  interference  modulation  function  allows  a  two 
parameter  search  using  data  from  several  narrowband  lines  in  a  manner  only 
dependent  on  the  underlying  source-receiver  geometry. 

The  simulated  annealing  algorithm  lends  itself  to  an  automatic  optimisation.  Being 
essentially  a  random  sampling  of  the  parameter  space,  it  is  not  an  efficient  technique 
however  it  is  robust.  We  find  that  such  an  algorithm  will  detect  the  correct  solution 
point  from  a  complex  cost  function  surface,  although  owing  to  the  need  to  evaluate 
many  cost  function  values,  it  can  be  slow.  This  is  the  most  straightforward  of  the  two 
methods,  and  by  choice  the  preferred  option.  Being  a  probabilistic  solution,  the 
solution  point  can  never  be  guaranteed  to  be  the  global  minimum,  and  the  results  from 
the  annealing  algorithm  always  should  be  verified  by  evaluation  of  the  modulation 
functions  for  each  of  the  narrowband  lines.  Any  significant  mismatch  would  indicate  a 
poor  solution. 

In  the  end,  the  decision  on  how  to  proceed  will  be  determined  by  the  quality  of  the 
data,  the  a  priori  knowledge  of  the  system  and  the  speed  of  computation. 


6.  References 

1)  P.R.  Lewis,  The  Effect  of  Ray  Curvature  on  Lloyd's  Mirror  Fringe  Estimates  and 
Source  Localisation,  DSTO  Technical  Report  DSTO-TR-0064, 1994. 

2)  SIMDOP16,  A  programme  to  produce  simulated  data  time  series,  D.R.A.  McMahon, 
unpublished. 

3)  AMADEUS,  Users'  Guide  V1.7,  1994,  D.R.A.  McMahon  &  G.K.  Schwarz,  DSTO 
General  Document  DSTO-GD-0073. 

4)  CASCADE,  A  Computer  Aided  Sonar  Classification  and  Data  Extraction  System,  Ebor 
Computing,  unpublished. 

5)  MATLAB  V4.2c,  1995,  The  Mathworks  Co. 

6)  The  1985  Baseline  Passive  RAYMODE  propagation  Loss  Prediction  programme 
(V8.0),  unpublished. 


21 


DSTO-TR-0583 


22 


DSTO-TR-0583 


Appendix  A 

Parameter  estimation  using  a  Simulated  Annealing  algorithm 

For  cost  function  surfaces  with  a  complex  structure  of  local  minima  and  valleys,  it  is 
non-trivial  to  produce  an  automated  optimisation  algorithm  which  is  capable  of 
locating  the  global  minimum  on  the  surface.  Standard  searching  algorithms  provide  a 
method  of  descending  a  local  valley,  generally  utilising  gradient  information,  and 
efficiently  locating  the  local  minimum  point,  and  are  not  suited  to  such  a  problem.  The 
simulated  annealing  algorithm  provides  a  way  of  sampling  the  parameter  space  in 
such  a  manner  that  it  will  converge  to  the  global  minimum  with  a  high  probability. 
The  method  requires  only  cost  function  value  evaluations  to  be  made,  and  as  such  is 
well  suited  to  the  optimisation  of  discrete  functions. 

The  algorithm  is  derived  by  analogy  with  the  adiabatic  cooling  of  a  thermodynamic 
system,  in  which  the  global  minimum  energy  state  is  achieved.  Adiabatic  cooling 
requires  that  the  temperature  is  reduced  sufficiently  slowly  that  the  system  remains  in 
thermal  equilibrium.  At  any  temperature,  the  particle  energies  are  distributed 
according  to  a  Boltzmann  distribution  and  there  is  always  some  probability  that  a 
particle  of  any  particular  energy  will  exist,  with  this  probability  diminishing  with  the 
reduction  in  temperature.  The  analogy  to  a  function  minimisation  proceeds  as  follows: 

•  The  system  energy  <=>  The  cost  function 

•  The  configuration  of  the  particles  o  The  set  of  possible  parameter  combinations 

•  The  temperature  <=>  The  searching  control  parameter 

The  basic  structure  of  the  algorithm  is  that  at  each  temperature,  a  random  sequence  of 
parameter  configurations  is  generated.  A  particular  configuration  is  accepted  if  it 
either  produces  a  smaller  than  current  cost  function  value,  or  with  a  Boltzmann 
distributed  probability  determined  by  the  temperature.  In  this  way  a  sampling  of  the 
parameter  space  is  performed,  and  after  some  number  of  configuration  steps,  the 
sampled  cost  function  should  converge  such  that  the  average  value  is  unchanging.  The 
best  point  found  during  this  sampling  then  serves  as  an  initial  point  for  the  next  cycle, 
which  occurs  at  a  lowered  temperature.  The  best  point  should  converge  to  the  optimal 
point  with  the  reduction  in  temperature.  Thus,  in  the  early  stages,  at  a  "high" 
temperature,  a  random  sampling  of  the  space  is  done.  As  the  temperature  is  reduced, 
the  best  point  will  settle  into  one  of  the  many  valleys,  although  it  is  always  possible  for 
this  point  to  jump  to  a  new  valley,  and  thus  allow  the  exploration  of  a  new  parameter 
region.  The  probability  of  doing  this,  however,  diminishes  with  the  temperature 
reduction. 

The  components  of  the  algorithm  which  need  to  be  specified  are  the  rules  for  the 
selection  of  a  parameter  configuration,  the  number  of  configurations  to  be  tested,  the 
method  of  reducing  the  temperature  at  each  cycle  and  the  termination  criteria.  Because 
die  algorithm  only  uses  the  cost  function  values,  and  samples  the  parameter  space  in  a 
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pseudo-random  manner,  it  will  generally  be  more  wasteful  of  function  evaluations 
than  one  guided  by  gradient  values,  however  this  also  makes  it  better  suited  to 
automated  optimisation.  The  details  of  the  algorithm  used  in  the  present  application 
are  given  below. 

As  observed  on  §3.1,  the  effect  of  the  reflection  coefficient  is  approximately  decoupled 
from  that  due  to  the  source  depth  and  range  at  CPA.  Therefore,  only  a  two  parameter 
search  is  considered,  and  the  cost  function  is  discretised  on  a  grid  constructed  by  die 
methods  of  §2.1.  An  initial  point  on  the  grid  is  selected  randomly. 

The  temperature  dependence  at  cycle  (i+1)  is  given  by 

Ti+1=aTi  with  0<a<l  (Al) 

The  value  of  a  will  determine  the  rate  at  which  the  system  is  cooled.  The  cooling  rate  is 
important  since  if  it  is  too  slow,  the  number  of  cycles  required  for  proper  convergence, 
and  hence  the  number  of  function  evaluations,  will  become  excessive.  Conversely,  if 
too  fast,  the  system  will  never  achieve  an  equilibrium  state  and  the  solution  point,  if 
any,  will  not  be  the  global  minimum.  The  appropriate  rate  of  cooling  will  depend  on 
the  scheme  used  to  select  the  parameter  configurations.  The  cooling  rate  can  be  larger 
if  the  potential  configurations  available  from  the  current  configuration  (the  "reach") 
cover  most  of  the  valid  parameter  space,  because  at  each  step  in  a  cycle  the  probability 
of  moving  out  of  the  current  region  to  a,  perhaps,  more  favourable  region  is  higher. 

It  is  desirable  to  bias  the  reach  toward  the  current  point,  since  then  more  of  the 
sampling  is  done  about  what  is  currently  the  best  estimate  of  the  global  minimum. 
Therefore,  within  a  cycle,  the  configuration  selection  is 

xi+i  =xi+Ax-sign(R[-l,l])R[0,l]N  with  Ax  =  P(xmax -xmin)  (A2) 

where  R[a,b]  is  a  uniformly  distributed  random  number  in  the  interval  [a,b]  and  Xi  is 
the  parameter  value  at  the  i*  step  of  the  cycle.  The  Ax  is  the  maximum  allowable  step. 
When  there  are  multiple  parameters,  each  parameter  is  sequentially  modified.  If  the 
new  parameter  value  exceeds  the  bounds  of  the  parameter  interval,  the  value  is 
reflected  back  into  the  interval  as 


x;.,— » 


IZx^-x^,  xi+1>Ximx 
^x^-x^,  Xitl<x„, 


(A3) 


The  non-uniform  random  distribution  of  the  step  sizes  allows  coverage  of  the  entire 
range  specified  by  Ax,  but  favours  configurations  near  to  the  current  position.  On  a 
discrete  grid,  the  value  of  (Xi+1)  is  rounded  to  match  the  nearest  grid  position.  Care 
needs  to  be  taken  if  the  discrete  grid  is  small,  as  the  power  N  in  equation  A2,  if  too 


24 


DSTO-TR-0583 


large,  will  cause  most  of  the  new  parameter  configurations  to  match  the  current 
configuration. 

A  parameter  configuration  is  accepted  if  either  the  resultant  cost  function  is  smaller 
than  the  current  cost  function,  or  if,  at  cycle  i,  the  change  in  the  cost  functions  AC 

satisfies 


/ 

R[0,1]  <  exp 

v 


AC' 

tJ 


for  AC  >  0 


(A4) 


The  initial  temperature  To  is  selected  to  be  larger  than  the  largest  expected  cost 
function  deviation. 


A  cycle  is  terminated  when  either  the  number  of  new  configurations  exceeds  a  preset 
level,  Ns,  or  the  number  of  attempts  at  finding  a  new  configuration  exceeds  a  preset 
level,  Na.  The  former  condition  carries  the  assumption  that  after  Ns  successful  new 
configurations,  the  system  has  settled  to  an  equilibrium  point.  The  latter  condition 
ensures  that  each  cycle  will  terminate  even  when  a  new  configuration  can  not  be 
found.  Appropriate  values  for  Ns  and  Na  will  depend  on  the  rate  at  which  the  system 
is  cooled.  At  the  termination  of  each  cycle  the  configuration  producing  the  lowest  cost 
function  is  chosen  as  the  initial  point  for  the  next  cycle  of  the  process,  continuing  at  a 
reduced  temperature. 

The  annealing  process  is  terminated  when  either  the  number  of  successful  new 
configurations  for  the  previous  cycle  is  zero,  the  same  best  configuration  is  repeatedly 
found  over  a  number  of  cycles.  Nr,  or  when  the  number  of  temperature  reductions 
exceeds  a  preset  level,  Nt-  The  first  condition  signifies  a  true  convergence  to  a  solution 
point;  similarly  with  the  second,  although  this  prevents  any  oscillation  about  the 
solution  point  over  several  cycles.  Again,  the  latter  condition  is  a  fail-safe  termination 
condition. 

The  values  chosen  for  the  present  application  are  Ns  =  15,  Na  =  60,  Nr  =  10  and 
Nt  =  40.  The  cooling  factor  is  a  =  0.65,  the  step  size  scale  factor  is  P  =  0.5,  and  the  step 
non-linear  random  number  distribution  is  generated  with  a  power  of  N  =  2.  Once 
suitable  choices  for  these  parameter  values  have  been  made,  the  annealing 
optimisation  is  particularly  straightforward.  The  difficulty  is  in  selecting  values  that 
will  ensure  an  adequate  sampling  of  the  parameter  space,  whilst  providing  good 
convergence  within  a  reasonable  amount  of  time. 

Examples  of  the  result  from  a  two  parameter  search  for  the  examples  of  §3.1  are  shown 
in  figures  14  and  15.  Figure  14  shows  the  sampled  cost  function  during  the  search,  and 
how  the  algorithm  samples  the  entire  parameter  space,  but  concentrates  on  the  region 
near  to  the  estimated  solution  point.  The  variation  in  the  cost  function  at  the  best 
estimate  of  the  solution  for  each  cycle,  and  examples  of  the  sampling  within  each  cycle 
are  shown  in  figure  15. 
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Appendix  B 


The  LINAL  analysis  interface 

To  facilitate  the  analysis  of  narrowband  data,  a  prototype  analysis  programme  has 
been  developed.  This  is  a  menu  driven  interface  written  using  the  MATLAB  [5]  data 
analysis  package  running  under  UNIX.  It  provides  all  of  the  analysis  features 
discussed  in  the  body  of  the  report,  and  was  used  in  analysing  the  simulation  data 
therein.  Most  of  the  functions  are  written  as  MATLAB  shell  scripts,  and  thus  are  easily 
modified.  To  improve  the  execution  speed,  the  cost  function  evaluation  routine  is  a 
compiled  C  module. 


COMMAND  STATUS 

SOUNDSPEED  PROFILE 
FILENAME 

LINE  DATA 
FILENAME 

BACKGROUND  DATA 
FILENAME 

WORKSPACE 

FILENAME 

LOSS  DATA 
FILENAME 

•  7  v  '  v  77  •  . 

.  •'  '  :• 

FILE 

OPTIONS 

CALCULATION 

OPTIONS 

DISPLA 

OPTION 

Y 

rs 

SOURCE 

PATH 

GRID 

CONTROL 

ANNEALING 

CONTROL 

MINIMISATION 

CONTROL 

GRID 

SOLUTION 

ANNEALING 

SOLUTION 

MINIMISATION 

SOLUTION 

The  menu  display  is  shown  in  figure  Bl.  The  information  is  divided  into  the  following 
sections: 


•  COMMAND  STATUS 

•  FILENAME 

•  FILE  OPTIONS 

•  CALCULATION  OPTIONS 

•  DISPLAY  OPTIONS 

•  SOURCE  PATH 

•  GRID  CONTROL 

•  ANNEALING  CONTROL 

.  MINIMISATION  CONTROL 

.  GRID  SOLUTION 

•  ANNEALING  SOLUTION 

•  MINIMISATION  SOLUTION 


the  status  of  the  current  command 
the  names  of  the  various  active  data  files 
the  file  open,  load  and  save  options 
the  various  data  calculation  options 
the  various  data  display  and  printing 
options 

the  description  of  the  current  source 
trajectory  solution 

the  parameter  grid  control  parameters 
the  simulated  annealing  optimisation 
control  parameters 

the  simplex  least  squares  optimisation 
control  parameters 
the  current  search  parameter  grid 
the  simulated  annealing  optimisation 
trajectory  solution 

the  simplex  least  squares  optimisation 
trajectory  solution 
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A  description  of  the  various  options  available  in  each  section  is  given  below. 

B.l  Filename 

Five  files  are  used  by  the  programme:  the  sound  speed  profile  description,  the 
narrowband  line  profile  data  file,  a  background  noise  spectrum  data  file,  a  workspace 
file  and  a  loss  data  file.  The  files  provide  for 

•  the  sound  speed  profile  description,  which  is  an  ASCII  text  file  containing  a  list 
of  depth,  sound  speed  pairs.  Each  record  lists  one  pair  of  space  separated  values 
giving  the  depth  (metres)  and  sound  speed  (metres/ second).  Blank  lines  and 
lines  with  a  %  (comment  lines)  as  the  first  character  are  ignored.  If  only  a  single 
data  pair  is  present,  it  is  interpreted  as  the  water  column  depth  and  the  constant 
sound  speed  for  this  column. 

•  the  line  profile  data,  which  is  a  CASCADE  [4]  formatted  text  file.  This  contains 
the  narrowband  line  Doppler  shifted  frequency,  signal  strength  and  time  values. 
It  is  produced  using  AMADEUS  [3]. 

•  the  complete  state  of  the  current  analysis  environment  to  be  saved  to  a 
MATLAB  binary  format  workspace  file.  This  allows  standard  configurations  to 
be  established,  or  for  a  continuation  of  an  analysis  session  at  a  later  date.  The 
currently  displayed  figures  are  not  saved. 

•  a  background  noise  level  frequency  spectrum  to  allow  background  noise 
subtraction.  The  file  is  an  ASCII  text  file  giving  the  power  level  in  dB/Hz.  Each 
line  contains  one  pair  of  space  separated  values  giving  the  frequency  (Hz)  and 
background  level  (dB/Hz).  Blank  lines  and  lines  with  a  %  (comment  lines)  as 
the  first  character  are  ignored. 

•  the  internally  generated  loss  data  to  be  saved  to  a  MATLAB  MAT  binary  file. 
This  may  be  of  use  when  the  required  loss  data  is  computationally  expensive  to 
produce. 

B.2  File  Options 

These  options  control  the  loading  and  saving  of  the  various  data  files.  The  options 
available  are  to 

•  load  the  sound  speed  profile 

•  load  the  narrowband  line  data 

•  load  a  background  spectrum 

•  load  previously  saved  loss  data 

•  load  a  previously  saved  workspace  description 

•  save  the  currently  evaluated  loss  data 

•  save  the  current  workspace  description 

•  exit  the  programme 

Loading  any  of  the  sound  speed,  line  data,  background,  workspace  or  loss  data  files 
will  cause  the  current  filename  entry  to  be  modified.  If  the  workspace  is  altered  from 
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that  defined  by  the  currently  loaded  file,  die  filename  is  marked  with  an  asterisk  to 
indicated  this.  When  the  line  data  file  is  loaded,  the  programme  will  prompt  for  the 
natural  frequency  of  each  of  the  narrowband  lines,  whilst  indicating  the  estimate  of 
this  based  on  the  doppler  shifted  frequency  values. 

B.3  Calculation  Options 

These  options  control  the  calculation  of  the  various  quantities,  and  the  operation  of  the 
optimisation  routines.  The  options  available  are  to 

•  generate  the  data  grid 

•  calculate  the  propagation  loss  data  on  the  defined  grid 

•  calculate  the  cost  function  on  the  defined  grid 

•  determine  an  optimal  solution  using  the  simulated  annealing  algorithm 

•  determine  the  optimal  solution  using  a  simplex  least  squares  minimisation 
algorithm 

•  set  the  data  smoothing/filtering  window  sizes  for  each  narrowband  line 

•  select  new  grid  limits  from  a  currently  displayed  cost  function  surface 

•  set  or  modify  the  background  levels  at  each  of  the  line  data  frequencies 

A  menu  button  is  used  to 

•  set  the  cost  function  weighting  mode 

Several  switches,  one  per  line  frequency,  to 

•  select  or  unselect  the  various  line  data  sets  from  the  optimisation 

The  sequence  for  performing  an  analysis  is  to  load  the  appropriate  sound  speed 
profile,  narrowband  line  data  and,  possibly,  a  background  noise  spectrum.  When  the 
line  data  file  is  loaded,  the  pushbuttons,  one  per  frequency,  will  be  enabled.  When 
pressed,  these  buttons  include/ exclude  the  associated  data  from  the  optimisation 
process. 

Before  any  calculation  can  proceed,  a  data  grid  must  be  defined  using  the  Data  Grid 
option.  The  size  of  the  grid  is  determined  by  the  parameters  in  the  Grid  Control 
control  panel.  The  size  of  the  grid  will  in  part  be  determined  by  the  highest  frequency 
selected,  and  increases  as  die  square  of  this  frequency.  All  calculations  are  based  on 
the  current  grid.  Changing  the  selected  frequencies  will  issue  a  warning  that  the 
displayed  grid,  loss  data  and  cost  function  data  are  no  longer  valid. 

The  Loss  Data  option  allows  the  loss  data  to  be  generated  for  the  current  grid  and  the 
current  frequency  selection.  Changing  the  frequency  selection  will  clear  the  current 
loss  data.  The  sound  speed  profile  determines  the  manner  of  the  propagation  loss 
calculation.  If  the  sound  speed  is  constant,  the  analytic  expressions  are  used, 
otherwise,  a  range  independent  propagation  loss  programme  is  utilised.  The 
propagation  loss  is  internally  generated  if  an  iso-speed  profile  is  used,  otherwise  via 
an  external  programme. 


29 


DSTO-TR-0583 


The  relevant  function  is  loss.m.  As  currently  setup,  this  calls  a  further  function 
raymode,  m  which  acts  as  an  interface  to  an  independent  propagation  programme 
RAYMODE  [6].  This  programme  utilises  the  UNIX  stdin  and  stdout  I/O  files  to  accept 
data  and  to  return  the  results  of  any  calculations.  The  UNIX  pipe  functionality  is  used 
within  raymode.m  to  send  a  formatted  input  deck  to  RAYMODE  and  to  read  the 
formatted  output.  A  UNIX  Awk  script  is  used  to  scan  the  programme  output  and 
extract  the  relevant  propagation  loss  data.  The  details  of  this  are  hidden  inside  the 
raymode.m  function.  This  function  and  die  calling  sequence  within  loss.m  is  specific  to 
the  actual  propagation  loss  programme  used.  Any  programme  that  uses  the  stdin  and 
stdout  for  data  1/ O  can  be  used,  although  the  function  to  call  it  and  extract  the  relevant 
data  will  need  to  be  tailored  to  the  programme.  The  calling  sequence  inside  loss.m  will 
also  need  to  be  altered. 

The  total  cost  function  surface  can  be  evaluated  for  the  current  frequency  selection  and 
the  current  grid  with  the  Gridded  Cost  Function  option.  The  loss  data  must  be 
defined  before  the  cost  function  is  evaluated.  Changing  the  frequency  selection  will 
clear  the  current  cost  function  data. 

The  simulated  annealing  solution,  using  the  parameter  values  currently  defined  in  the 
Annealing  Control  panel,  is  produced  using  the  Annealed  Solution  option.  Only  the 
source  depth  and  the  range  at  CPA  are  determined.  No  uncertainty  estimates  are 
produced  All  trajectory  and  receiver  parameters  are  obtained  from  those  currently 
defined  in  the  Source  Path  control  panel. 

The  least  squares  optimisation,  using  the  parameters  currently  defined  in  the 
Minimisation  Control  panel,  is  produced  using  the  Fit  Solution  option.  The  source 
depth,  range  at  CPA  and  the  reflection  coefficient  are  determined,  as  specified.  All 
trajectory  and  receiver  parameters  are  obtained  from  those  currently  defined  in  the 
Source  Path  control  panel.  A  simplex  algorithm,  which  uses  only  the  function  values 
at  the  simplex  vertices,  is  used.  No  constraints  are  imposed  on  the  parameter  values. 
The  MATLAB /mms.m  function  implements  this  function  minimisation  routine. 

The  definition  of  the  cost  function  requires  an  estimate  of  the  unmodulated  line 
profiles.  These  are  determined  as  discussed  in  §2.3,  and  require  an  estimate  of  the 
smoothing  filter  window.  The  size  of  this  will  depend  on  the  frequency  of  the 
narrowband  line  data.  The  estimates  are  currently  obtained  from  user  defined  values. 
The  Smoothing  Window  option  is  used  for  this.  On  selection,  each  line  data  set  will 
be  displayed  in  sequence.  The  positions  of  the  interference  minima  must  be  marked  in 
using  the  cursor  (left  button  to  select,  middle  button  to  remove,  right  button  to  end).  A 
window  size  is  then  determined  from  an  average  of  the  marked  minima  separations 
about  the  time  of  CPA.  It  is  therefore  important  that,  where  possible,  adjacent  minima 
are  marked.  For  reference,  the  current  value  of  the  time  of  CPA  defined  in  the  Source 
Path  control  panel  will  be  shown.  The  estimated  window  widths  at  each  frequency  are 
listed. 

The  Select  on  Grid  option  is  provided  to  simplify  the  selection  of  a  region  of  the 
parameter  space  of  interest,  and  can  be  selected  whenever  a  cost  function  surface  is 
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displayed.  A  sub-region  of  the  current  parameter  space  (the  source  depth  and  die  CP  A 
range)  is  selected  using  the  cursor.  The  appropriate  limits  are  then  transferred  to  the 
Grid  Control  panel.  To  operate  the  cursor,  firstly  a  figure  displaying  a  cost  function 
surface  is  selected  (left  button),  followed  by  the  sub-region  definition  (left  button).  A 
rubberband  box  shows  the  selected  region. 

The  background  noise  levels  can  be  optionally  included  via  the  Load  NB  Line  Data 
option  of  the  File  Option  panel.  If  no  file  is  loaded,  a  zero  level  background  is 
assumed.  Alternatively,  the  Set  Background  Levels  option  can  be  used.  This  will 
modify  the  current  levels,  be  they  externally  defined,  or  the  default  levels.  The  values 
are  in  dB/Hz,  and  are  specified  at  the  defined  line  data  frequency  values. 

The  cost  function  is  evaluated  from  the  set  of  line  data  signal  levels.  The  contribution 
of  each  datum  can  be  optionally  weighted.  The  No  Weighting/Smoothing  menu 
button  determines  the  current  state.  The  data  weighting,  if  applied,  follows  the 
definition  of  §2.3.  By  preference,  the  data  should  be  weighted. 

Switches  control  the  selection  of  the  data  to  be  included  in  any  optimisation.  Changing 
the  data  selection  will  render  any  previously  calculated  propagation  loss  data, 
parameter  grid  and  cost  function  values  invalid. 

B.4  Display  Options 

A  number  of  display  and  hard  copy  options  are  provided.  These  are  collected  in  the 
Display  Options  panel,  and  comprise  options  to 

•  plot  the  current  sound  speed  profile 

•  plot  all  the  currently  loaded  narrowband  line  data 

•  plot  all  the  currently  loaded  line  data  together  with  the  current  fit 

•  plot  the  fit  residual  based  on  the  current  line  data  and  fit 

•  plot  the  smoothed  current  line  data 

•  plot  the  ratio  of  the  current  line  data  to  the  current  smoothed  line  data 

•  plot  the  current  propagation  loss  data 

•  plot  the  current  cost  function  generated  with  the  Gridded  Cost  Function 
option 

•  plot  the  current  cost  function  resulting  from  the  current  simulated  annealing 
optimisation 

•  print  the  figure  as  selected 

•  print  the  current  set  of  figures 

•  clear  the  figure  as  selected 

•  clear  the  current  set  of  figures 

The  number  of  plots  shown  for  any  of  the  plotting  options  depends  on  the  data  to  be 
displayed.  The  current  set  of  figures  is  consists  of  the  complete  set  of  plots  created  by 
the  most  recent  plotting  command.  This  set  can  be  printed,  with  one  plot  per  page,  or 
deleted.  In  this  latter  case,  the  current  set  of  figures  reverts  to  the  previously  most 
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recently  defined  set.  As  well,  an  individual  plot  can  be  printed  or  deleted.  When 
chosen,  a  prompt  for  the  figure  to  be  printed  or  deleted  is  given,  with  the  selection 
being  made  with  the  cursor. 

The  Sound  Speed  Profile  plot  consists  of  a  single  figure  showing  the  current  sound 
speed  profile. 

The  NB  Line  Data  plot  consists  of  multiple  figures,  each  showing  six  line  data  sets. 
All  the  currently  defined  data,  irrespective  of  the  status  of  the  current  data  selection, 
are  displayed.  Similarly,  the  NB  Line  Data  &  Fit  plot  consists  of  multiple  figures, 
each  with  six  line  data  sets,  together  with  the  estimated  profiles  based  on  the  current 
trajectory  parameters.  The  loss  data  necessary  for  these  profiles  are  calculated  when 
this  option  is  chosen,  with  the  manner  of  calculation  based  on  the  current  sound  speed 
profile.  The  data  in  both  cases  are  plotted  in  absolute  units  (ie.  not  dB). 

The  NB  Fit  Residual  plot  consists  of  multiple  figures,  each  showing  six  line  data  sets. 
Here,  the  difference  between  the  line  data  and  the  calculated  profiles  is  shown.  As 
above,  the  loss  data,  using  the  current  trajectory  parameters,  are  calculated  when  this 
option  is  chosen.  The  data  are  plotted  in  dB. 

The  NB  Smoothed  Data  plot  again  consists  of  multiple  figures,  with  six  line  data  sets 
shown  per  plot.  In  this  case  the  data  before  and  after  being  passed  through  the 
smoothing  filter  are  displayed.  Similarly,  the  ratio  of  the  line  data  to  the  smoothed 
version  of  the  same  data  is  shown  with  the  NB  Data/Smoothed  Data  option.  These 
latter  plots  are  in  linear  units,  and  if  the  data  smoothing  is  adequate,  should  show  as 
an  oscillatory  function  of  roughly  constant  amplitude.  The  size  of  the  oscillations 
depends  on  the  reflection  coefficient,  and  as  a  guide  to  a  suitable  value,  the  lower  and 
upper  bounds  of  the  oscillation,  evaluated  for  an  iso-speed  profile  (Cs  =  1500  m/s)  and 
with  the  reflection  coefficients  given  in  the  Source  Path  panel,  are  plotted. 

If  the  current  propagation  loss  data,  calculated  using  the  Loss  Data  option  in  the 
Calculation  Options  panel,  exists,  it  is  plotted  using  the  Loss  Data  option.  One  figure 
is  shown  for  each  active  line  data  set.  The  loss  data  are  plotted  as  a  function  of  the 
source  depth  and  the  horizontal  range  from  the  source  for  a  fixed  receiver  depth  and 
source  frequency. 

The  Gridded  Cost  Function  option  plots  the  total  cost  function  calculated  on  the 
parameter  grid  using  the  Gridded  Cost  Function  option  in  the  Calculation  Options 
panel.  Only  the  currently  active  line  data  sets  are  used  in  the  cost  function  evaluation. 
A  single  figure  is  produced.  If  a  simulated  annealing  solution  has  been  obtained,  a  by¬ 
product  is  the  total  cost  function  evaluated  at  the  sampled  parameter  points.  This  can 
be  used  to  show  how  the  parameter  space  has  been  covered,  or  to  define  a  smaller 
region  of  interest.  The  Annealed  Cost  Function  option  plots  this  function,  if  it  exists. 
A  single  figure  is  produced.  The  cost  function  for  both  these  options  is  plotted  as  a 
function  of  the  source  depth  and  the  horizontal  range  at  CPA. 
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B.5  Source  Path 

The  current  source  straight  line  trajectory  is  specified  using  this  panel,  and  all 
calculations  are  based  on  the  current  values  of  the  parameters.  These  parameters  are 
the 

•  horizontal  CPA  range  (metres) 

•  slant  CPA  range  (metres) 

•  time  of  CPA  (seconds) 

•  source  speed  (metres/ second) 

•  source  depth  (metres) 

•  receiver  depth  (metres) 

•  surface  reflection  coefficient 

•  bottom  reflection  coefficient 

•  smoothing  factor 

The  quantities  are  self-evident.  The  CPA  range  can  be  entered  as  either  the  horizontal 
or  the  slant  range.  The  associated  value  will  be  calculated  from  the  current  range  and 
depth  parameter  values.  The  reflection  coefficients  are  plane  wave  frequency  and 
angle  independent  values.  A  reflection  coefficient  of  zero  will  disable  reflection  from 
the  associated  surface. 

The  current  values  of  these  quantities  are  used  in  producing  the  modelled  line  profiles 
(ie.  the  loss  data),  in  the  determination  of  the  parameter  grid  size,  for  defining  the 
initial  parameter  values  for  the  simulated  annealing  minimisation,  and  for  the  constant 
parameter  values  in  the  simplex  least  squares  minimisation. 

The  values  for  the  CPA  range  and  the  source  depth  can  be  read  from  the  current  cursor 
position  on  a  cost  function  plot  using  the  Cursor  Selection  button.  Firstly,  a  figure 
displaying  a  cost  function  surface  is  selected  with  the  cursor  (left  button).  The  cursor  is 
then  position  on  the  surface  and  a  point  selected  (left  button).  This  point  will  be 
indicated  with  a  cross  hair.  The  values  are  entered  in  the  Source  Path  panel.  This 
allows  primarily  an  estimate  of  the  position  of  a  local  minimum  to  be  read  and  the 
resultant  match  to  the  data  to  be  tested.  This  can  be  used  to  assess  the  quality  of  the 
local  minimum  of  the  cost  function  as  the  best  solution. 

B.6  Grid  Control 

All  calculations  are  based  on  a  discrete  parameter  grid.  This  panel  provides  the  control 
parameters  for  the  grid  determination.  Three  discrete  sets  are  calculated,  for  the  source 
depth  and  the  horizontal  range  at  CPA,  and  for  the  horizontal  ranges  associated  with 
the  propagation  loss  calculations.  The  parameters  are  the 

•  lower  and  upper  limits  on  the  horizontal  CPA  range  (metres) 

•  lower  and  upper  limits  on  the  source  depth  (metres) 

•  step  size  scale  factor  for  the  range  step  in  the  loss  calculations 

•  step  size  scale  factor  for  the  horizontal  CPA  range 

•  step  size  scale  factor  for  the  source  depth 
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As  well,  three  switches  are  provided  to  restrict  the  step  sizes  to 

•  a  constant  loss  range  step 

•  a  constant  CPA  range  step 

•  a  constant  source  depth  step 

The  lower  and  upper  limits  on  the  CPA  range  and  the  source  depth  define  the  region 
of  parameter  space  of  interest.  The  range  limits  for  the  propagation  loss  data  are 
determined  from  the  CPA  range,  the  time  of  CPA  and  the  time  limits  of  the  line  data. 
The  determination  of  the  grid  is  defined  in  equations  7  and  9,  which  give  the 
maximum  step  sizes  consistent  with  the  non-aliasing  of  the  data  within  an  iso-speed 
model. 

The  scale  factors  are  needed  for  three  reasons.  To  produce  the  estimates  of  the  line 
profile  data,  the  loss  data  is  linearly  interpolated  at  the  ranges  defined  by  the  times 
and  the  current  path  parameter  values.  The  quality  of  this  interpolation  will  improve 
as  the  step  size  is  reduced.  Secondly,  the  necessary  step  sizes  are  expected  to  be  smaller 
when  a  non-constant  sound  speed  profile  is  used  (eg.  as  shown  by  the  lines  of  constant 
phase  on  a  lofargram)  [1].  In  this  case,  it  may  be  necessary  to  adjust  the  step  size 
empirically  in  order  that  aliasing  of  the  loss  data  does  not  occur.  Lastly,  there  may  be  a 
requirement  for  a  detailed  exploration  of  part  of  the  cost  function  surface. 

The  step  sizes  can  be  held  constant  at  the  smallest  calculated  values.  The  constant  step 
switches  select  this  option.  In  particular,  depending  on  the  programme  being  used  the 
propagation  loss  data  may  need  to  be  evaluated  at  equidistant  intervals;  thus  the  need 
for  this  functionality. 

The  current  grid  size  is  displayed  in  the  Grid  Solution  panel.  This  is  updated  using 
the  current  parameters  in  the  Grid  Control  panel  along  with  the  maximum  active  line 
frequency.  The  size  of  the  grid  and  the  minimum  and  maximum  step  sizes  for  each 
relevant  parameter  are  displayed. 

B.7  Annealing  Control 

This  panel  maintains  the  control  values  for  the  simulated  annealing  minimisation.  The 
parameters  are  described  in  appendix  A.  The  inputs  are  the 

•  maximum  number  of  iterations  per  variable  at  each  cycle  (NA) 

•  maximum  number  of  successful  configurations  found  per  variable  at  each  cycle 
(NS) 

•  maximum  number  of  cycles  in  the  search  (Nr) 

•  cooling  factor  (a) 

•  step  scale  factor  (3) 

•  step  size  power  (N) 
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Two  switches  are  provided  to 

•  list  the  information  on  the  search  as  is  proceeds 

•  plot  the  evolution  of  the  search 

The  parameters  correspond  to  those  described  in  appendix  A,  as  indicated  by  the 
quantities  in  the  parentheses. 

The  minimisation  information  listed  when  the  List  Fit  Data  option  is  selected 
includes  the 

•  iteration  number 

•  current  temperature 

•  current  CPA  range  and  source  depth 

•  number  of  successful  steps  in  the  current  cycle 

•  number  of  cost  function  evaluations  made 

At  the  end  of  the  search,  the  information  includes  the 

•  final  temperature 

•  final  CPA  range  and  source  depth 

•  number  of  cost  function  evaluations  made 

•  %  hit  rate  (As  the  cost  function  is  evaluated,  it  is  saved.  The  number  of  cost 
function  values  that  have  been  obtained  from  this  saved  list  is  expressed  as  a 
percentage  of  the  total  number  of  cost  function  evaluations  made) 

•  %  coverage  of  the  grid 

•  %  coverage  of  the  CPA  range  space 

•  %  coverage  of  the  source  depth  space 

Ideally,  the  %  coverage  in  the  CPA  range  and  source  depth  should  be  close  to  100%, 
whilst  for  a  large  grid  the  %  coverage  of  the  grid  should  be  small.  When  the  Display 
Fit  Data  option  is  selected,  the  information  is  shown  on  two  plots.  One  shows  the  cost 
function  for  each  successful  configuration  in  the  current  step.  The  initial  value  is 
marked  by  a  horizontal  line.  The  average  value  of  the  cost  function  is  overlaid  on  this. 
The  other  plot  shows  the  variation  in  the  best  cost  function  value  at  each  cycle,  and 
thus  how  well  the  system  in  converging. 

The  solution  point  is  displayed  in  the  Annealing  Solution  panel.  The  two  values  thus 
produced  can  be  transferred  to  the  Source  Path  panel  and  used  in  the  description  of 
the  current  trajectory  with  the  Transfer  Solution  option.  No  uncertainty  estimates  are 
produced.  The  simplex  least  squares  method  can  be  used  to  generate  these  by  using 
the  simulated  annealing  solution  as  the  initial  point  to  this  routine. 
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B.8  Minimisation  Control 

The  simplex  least  squares  optimisation  is  controlled  from  this  panel.  Four  parameters 
can  be  selectively  varied.  The  initial  values  for  these  must  be  set.  Three  optimisation 
control  values  are  required.  The  inputs  are  the 

•  horizontal  CPA  range  (metres) 

•  source  depth  (metres) 

•  surface  layer  reflection  coefficient 

•  bottom  layer  reflection  coefficient 

•  maximum  number  of  search  iterations 

•  approximate  fractional  tolerance  required  for  each  parameter 

•  required  absolute  cost  function  tolerance 

One  switch  is  provided  to 

•  list  the  information  on  the  search  as  is  proceeds 

A  switch  is  associated  with  each  of  the  four  parameters.  This  selects  or  unselects  the 
associated  parameter  from  being  varied  during  the  optimisation.  The  initial  values  for 
all  four  parameters  are  those  given  on  this  panel. 

The  optimisation  control  values  default  to  a  limit  of  200  iterations,  a  1%  fractional 
tolerance  in  each  parameter  and  an  absolute  tolerance  on  the  cost  function  of  0.0001. 
The  success  of  the  parameter  search  will,  in  particular,  depend  on  these  latter  two 
values.  A  successful  search  is  achieved  when  both  tolerance  requirements  have  been 
met.  The  appropriate  values  in  each  case  will  depend  on  the  structure  and  smoothness 
of  the  cost  function  surface;  viz.  the  absolute  magnitude  of  the  cost  function  values, 
and  its  variation  with  the  parameter  values. 

The  minimisation  information  listed  when  the  List  Fit  Data  option  is  selected 
includes  the  diagnostic  information  from  the  simplex  search,  being  primarily  the 
current  simplex  value. 

The  solution  point  is  displayed  in  the  Minimisation  Solution  panel.  Uncertainty 
estimates,  based  on  the  linearised  least  squares  model  outlined  in  §4,  are  produced  and 
displayed  with  the  solution  point. 

The  four  values  thus  produced  can  be  transferred  to  the  Source  Path  panel  and  used 
in  the  description  of  the  current  trajectory  with  the  Transfer  Solution  option. 
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Figure  1  Surfaces  of  the  source  depth  and  CPA  range  step  sizes  as  a  function  of  the  source 
depth  and  CPA  range.  The  values  are  obtained  from  equation  7,  and  are  for  a 
receiver  depth  of  100  m  and  a  frequency  of  150  Hz.  The  upper  figure  shows  the 
absolute  value  of  the  step  size  in  the  CPA  range,  and  the  lower  figure  the 
absolute  value  of  the  step  size  in  the  source  depth.  The  heavy  curves  on  these 
surfaces  mark  the  line  of  minimum  value,  as  required  by  the  data  grid. 
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Figure  2a  Examples  of  the  effect  of  smoothing  the  narroiuband  Lloyd's  mirror  interference 

pattern.  The  simulated  data  is  for  a  dipole  beam  pattern.  This  figure  shows  the 
effect  on  the  line  profile,  where  the  red  and  purple  lines  show  the  modelled  and 
estimated  beam  patterns  respectively.  The  50  Hz  data  ivith  only  a  single  node 
pair  is  poorly  reproduced,  ivhilst  the  150  Hz  data  with  tiuo  node  pairs  is 
marginal.  The  higher  frequency  data  are  well  reproduced. 
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Figure  3a 


The  effect  of  data  smoothing  ivith  a  change  in  the  phase  value  behveen  the 
direct  and  reflected  waves.  The  surface  reflection  coefficient  is  0.25.  The  purple 
curve  shoxvs  the  modelled  line  profile.  The  upper  plot  corresponds  to  a  phase 
difference  at  the  CPA  time  of  n,  the  middle  to  0.5 tv  and  the  lower  to  0.  The  log 
of  the  data  has  been  smoothed  ivith  a  window  four  times  the  ividth  of  the 
estimated  interference  fringe  ividth. 
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Figure  3b  The  effect  of  data  smoothing  ivith  a  change  in  the  phase  value  betxoeen  the 
direct  and  reflected  waves.  The  surface  reflection  coefficient  is  0.8.  The  purple 
curve  show s  the  modelled  line  profile.  The  upper  plot  corresponds  to  a  phase 
difference  at  the  CPA  time  of  n,  the  middle  to  0.5 n  and  the  loioer  to  0.  The  log 
of  the  data  has  been  smoothed  with  a  window  four  times  the  width  of  the 
estimated  interference  fringe  width. 
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Figure  4  The  gridded  least  squares  cost  surface  about  tlte  knoxvn  solution  point.  The  cost 
function  is  defined  in  equation  11.  The  modelled  source  has  an  isotropic  beam 
pattern.  A  single  narroxuband  line  of  frequency  of  150  Hz  was  used  and  no  data 
xoeighting  has  been  included.  The  geometry  is  described  in  the  text. 
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The  gridded  least  squares  cost  surface  about  the  known  solution  point.  The  cost 
function  is  defined  in  equation  11.  The  modelled  source  has  an  isotropic  beam 
pattern.  Six  narrowband  lines  of  frequencies  of 50, 150,  300,  400,  550  and  750 
Hz  were  used.  The  cost  function  data  has  not  been  weighted.  The  geometnj  is 
the  same  as  for  figure  4,  and  is  described  in  the  text. 


Figure  6a  The  modulation  functions  for  the  six  narrowband  lines  of  figure  5.  This  figure 
is  for  the  correct  surface  reflection  coefficient  of  0.5.  The  purple  curve  shoivs  the 
modelled  modulation  function.  The  red  bands  define  tlte  modelled  variation 
limits  on  the  modulation  function.  The  data  for  the  50  Hz  line  is  in  particular 
badly  represented. 
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The  modulation  functions  for  the  six  narrowband  tines  of  figure  5.  This  figure 
is  for  a  reduced  reflection  coefficient  of  0.3.  The  purple  curve  shows  the 
modelled  modulation  function.  The  red  bands  define  the  modelled  variation 
limits  on  the  modulation  function.  The  data  for  the  50  Hz  tine  is  in  particular 
badly  represented. 
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Figure  7  The  relationship  between  the  narrowband  line  frequencies  used  with  the  present 
sinndations  and  the  broadband  Lloyd's  mirror  interference  pattern.  The 
narrowband  frequencies  have  been  selected  so  that  varying  numbers  of 
interference  node  pairs  and  phases  at  the  CPA  time  are  associated  with  each 
narrowband  line.  The  heavy  vertical  lines  in  this  figure  represent  the 
narrowband  lines  and  the  light  solid  curves  the  contours  of  the  interference 
minima. 
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Figure  9a  The  effect  of  smoothing  on  the  line  profiles  for  the  six  narrowband  lines.  Tltis 
figure  is  for  the  symmetric  dipole  beam  pattern.  The  red  curve  shows  the 
modelled  line  profile,  ivhilst  the  purple  curve  shows  the  profile  estimated  from 
the  smoothing  filter.  The  50  Hz  beam  pattern  is  poorly  reproduced.  The  quality 
of  the  match  improves  as  the  number  of  node  pairs  in  the  line  profile  is 
increased.  In  each  case  the  smoothing  ivindoio  ivas  four  times  the  estimated 
interference  fringe  spacing. 
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Figure  9b  The  effect  of  smoothing  on  the  line  profiles  for  the  six  narrowband  lines.  This 
figure  is  for  the  asymmetric  dipole  beam  pattern.  The  red  curve  shoxvs  the 
modelled  line  profile,  whilst  the  purple  curve  shoxvs  the  profile  estimated  from 
the  smoothing  filter.  The  50  Hz  beam  pattern  is  poorly  reproduced.  The  quality 
of  the  match  improves  as  the  number  of  node  pairs  in  tlie  line  profile  is 
increased.  In  each  case  the  smoothing  windoxv  was  four  times  the  estimated 
interference  fringe  spacing. 
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Figure  10a  The  modulation  functions  deduced  from  the  symmetric  dipole  pattern  data 
shown  in  figure  9a.  The  patterns  for  this  and  asymmetric  dipole  patterns  in 
figure  10b  are  similar,  shoioing  that  the  effect  of  the  beam  pattern  has  been 
removed.  The  purple  curve  shows  the  modelled  modulation  function. 
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Figure  10b  The  modulation  functions  deduced  from  the  asymmetric  dipole  pattern  data 
"  shown  in  figure  9b.  The  patterns  for  this  and  symmetric  dipole  patterns  in 

figure  10a  are  similar,  showing  that  the  effect  of  tlxe  beam  pattern  has  been 
removed.  The  purple  curve  shoxvs  the  modelled  modulation  function. 
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Figure  11a  The  least  squares  cost  function  surfaces  for  a  single  narroivband  line  of  150  Hz. 

The  data  is  for  the  symmetric  dipole  beam  pattern.  The  upper  figure  does  not 
include  any  data  iveighting.  The  loiuer  figure  includes  data  iveighting.  The 
iveighted  cost  function  surface  shows  a  smoother  variation  and  the  minimum 
point  matches  the  correct  solution  minimum.  The  umveighted  surface  is  more 
complex,  and  in  this  example  wrongly  identifies  the  solution  minimum. 
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Figure  lib  The  least  squares  cost  function  surfaces  for  a  single  narrowband  line  of  150  Hz. 

The  data  is  for  the  asymmetric  dipole  beam  pattern.  The  upper  figure  does  not 
include  any  data  iveighting.  The  loiver  figure  includes  data  weighting.  These 
figures  are  similar  to  those  in  figure  11a. 
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Figure  12  The  result  of  a  global  iveighted  least  squares  minimisation  of  the  fit  to  the 
symmetric  dipole  data.  The  six  narrowband  lines  at  frequencies  of 50,  150,  300, 
400,  550  and  750  Hz  were  included  in  tlte  optimisation.  A  direct  simplex 
searching  algorithm  xoas  used.  The  initial  point  was  taken  from  the  minimum 
identified  in  figure  11a.  The  estimated  solution  point  is  at  (Ro,  Sd,  rj)  = 
(476.2,  49.9,  0.457),  in  comparison  with  the  true  point  at  (477,  50,  0.5). 
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Figure  13  The  estimated,  beam  patterns  for  tire  six  data  sets  of  figure  12.  The  pattern 
should  he  identical  in  shape  in  each  case.  The  red  curve  shoivs  the  true  beam 
pattern.  These  have  been  normalised  to  the  data  since  no  calibration  of  the  data 
was  included.  Good  agreement  in  the  true  and  extracted  line  shapes  is  obtained 
in  each  case. 
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Figure  14  The  sampled  weighted  least  squares  cost  function  as  produced  using  the 
simulated  annealing  algorithm.  The  correct  minimum  point  is  identified.  The 
statistics  on  this  optimisation  are  given  in  the  text.  The  standard  annealing 
parameters  described  in  Appendix  A  were  used.  Approximately  2%  of  the 
available  grid  has  been  evaluated.  This  figure  shows  how  the  entire  range  lias 
been  sampled ,  but  that  the  region  closest  to  the  solution  point  has  received  the 
more  detailed  inspection. 
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Figure  15  A  snapshot  of  the  simulated  annealing  search.  This  figure  shows  the  progress 
made  at  the  eighth  cycle  of  the  search.  The  left  figure  shows  the  sample  cost 
function  values  for  the  current  cycle,  with  the  horizontal  line  indicating  the 
current  best  cost  function  value.  A  better  solution  point  has  been  identified. 
The  right  curve  shows  the  progress  in  the  best  cost  function  value  with  cycle 
number. 
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Figure  16  Correlations  between  the  various  trajectory  parameters  obtained  from  a 
variational  study.  The  left  column  shows  the  variation  of  the  deduced  unknown 
parameters  with  the  known  source  speed.  The  right  column  shows  the  variation 
with  the  known  receiver  depth.  The  data  are  normalised  to  the  true  values  given 
in  the  text.  The  solid  curves  indicate  the  expected  correlation  based  on  equation 
17.  The  source  depth  and  CPA  range  appear  to  be  strongly  dependent  on  the 
source  speed.  The  reflection  coefficient  is  independent  of  the  geometric 
parameters. 
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The  variation  in  the  estimated  solution  point  with  changes  in  the  known 
receiver  depth  and  source  speed.  The  data  are  presented  relative  to  the  true 
values.  The  solid  curve  shows  the  dependence  expected  from  equation  17.  The 
large  change  in  the  solution  point  matches  the  large  change  imposed  on  the 
source  speed  (±15%). 
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Figure  18  Estimated  error  bounds  on  the  solution  point  for  the  symmetric  dipole  beam 
pattern.  The  data  for  the  single  narroiuband  line  at  150  Hz  were  included.  The 
figure  shows  the  iveighted  cost  function.  The  single  contour  represents  the  95% 
confidence  level  based  on  equation  18.  Tlie  rectangular  box  shozos  the 
independent  uncertainty  95%  confidence  bounds,  based  on  the  linear  least 
squares  approximation  given  in  equations  19-22. 
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Figure  B1 


The  FINAL  MATLAB  demonstration  interface.  The  various  control  buttons 
and  options  are  described  in  Appendix  B. 
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